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QUANTITATIVE  VALIDATION  OF  A  MODEL  OF  CHLORINATED 
ETHENE  NATURAL  ATTENUATION 


1.0  Introduction 

1.1  Motivation 

The  contamination  of  groundwater  by  chlorinated  solvents  is  a  well-recognized 
problem  (Mackay  and  Cherry,  1989;  Adamson  and  Parkin,  2000).  Two 
chemicals  in  particular  pose  a  significant  threat  to  human  health  and  the 
environment:  tetrachloroethylene  (PCE),  and  trichloroethylene  (TCE).  This 
research  will  focus  on  the  treatment  of  aquifers  containing  PCE  and  TCE. 

Both  PCE  and  TCE  are  suspected  carcinogens  (Sturchio  et  al.,  1998),  making 
groundwater  contaminated  with  these  compounds  a  health  hazard.  People 
exposed  to  TCE  have  reported  an  excess  number  of  adverse  health  effects,  to 
include  liver  and  kidney  disease,  diabetes,  and  stroke  (ATSDR,  1996). 
Accordingly,  the  United  States  Environmental  Protection  Agency  (USEPA)  has 
set  the  Maximum  Contaminant  Level  (MCL)  for  both  PCE  and  TCE  in  drinking 
water  at  five  parts  per  billion  (CFR,  2000b)  and  the  Maximum  Contaminant  Level 
Goal  (MCLG)  at  zero  parts  per  billion  (CFR,  2000a). 

The  health  risk  from  PCE  and  TCE  is  based  on  exposure  as  well  as  toxicity. 

PCE  and  TCE  are  common  industrial  solvents,  used  mostly  to  degrease  metal 
and  to  produce  inks  and  paints.  After  years  of  intentional  and  unintentional 
releases,  these  solvents  have  infiltrated  into  the  ground  and  contaminated 
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underlying  aquifers.  The  extent  of  contamination  across  the  country  is  such  that 
the  EPA  has  declared  TCE  and  its  degradation  products,  cis-1 ,2-dichloroethylene 
(DCE)  and  vinyl  chloride  (VC),  to  be  priority  pollutants  (Bloom  et  al.,  2000).  TCE 
is  the  most  frequently  detected  groundwater  contaminant  at  hazardous  waste 
sites  in  the  United  States  (Bloom  et  al.,  2000).  Of  the  nine  Marine  Corps  facilities 
on  the  National  Priorities  List,  eight  of  them  have  PCE  or  TCE  contamination  of 
soil  or  groundwater  (USEPA,  2000). 

Remediating  a  contaminated  aquifer  is  not  a  simple  task.  Both  PCE  and  TCE 
are  classified  as  dense  non-aqueous  phase  liquids  (DNAPLs).  DNAPLs  are 
denser  than  water,  and  tend  to  sink  after  reaching  the  water  table.  While  doing 
so,  some  of  the  DNAPL  partitions  into  the  aqueous  phase  and  is  carried  off  by 
groundwater  flow.  Due  to  the  relatively  low  solubility  of  PCE  and  TCE  in  water,  a 
DNAPL  source  can  persist  for  decades  (Mackay  and  Cherry,  1989).  While  some 
experiments  have  had  preliminary  success  (Ho  et  al.,  1999),  there  is  no 
conventional  technology  that  effectively  removes  DNAPL  sources  (NRC,  1999). 
This  research  will  focus  on  containment  and  removal  of  contaminants  from  the 
aqueous  phase.  At  present,  there  are  three  strategies  being  implemented  to  deal 
with  solvent-laden  groundwater:  conventional  pump  and  treat,  passive  barriers, 
and  natural  attenuation. 

Conventional  pump  and  treat  is  by  far  the  most  common  approach  to  the 
containment  and  removal  of  TCE  and  PCE.  In  this  strategy,  water  is  pumped  to 
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the  surface  where  contaminants  are  removed  using  an  engineered  treatment 
process.  The  treated  water  is  then  either  returned  to  the  subsurface  to  prevent 
aquifer  dewatering  and  to  sustain  beneficial  hydraulic  gradients  or  discharged  to 
surface  water.  The  aboveground  contaminant  treatment  technologies  are  well 
understood,  as  is  the  hydraulic  containment  achievable  through  the  use  of 
pumping  wells.  While  pump  and  treat  poses  an  increased  risk  to  receptors  as 
contaminants  are  pumped  to  the  surface,  the  largest  drawback  of  this  strategy  is 
its  high  cost.  Due  to  the  fact  that  pump  and  treat  systems  must  operate  for  many 
years  to  contain  a  contaminant  plume,  lifecycle  costs  for  these  systems  are  quite 
high.  These  high  costs,  as  well  as  the  increased  risk  to  potential  receptors,  have 
prompted  the  development  of  alternative  strategies  such  as  passive  barriers  and 
natural  attenuation  to  manage  contaminated  groundwater. 

Passive  barriers  have  received  increased  attention  as  an  effective,  low-cost 
containment  strategy.  By  constructing  a  trench  filled  with  zero-valent  iron  filings, 
a  barrier  is  formed  that  contaminated  groundwater  flows  through.  As  PCE  and 
TCE  are  transported  through  the  barrier,  they  are  reductively  dechlorinated  in  the 
trench.  While  this  technology  is  fairly  well  understood  and  cost  effective,  it  has 
its  limitations.  The  system  is  only  appropriate  for  relatively  shallow  aquifers,  as 
the  technology  is  limited  to  the  depth  a  trench  can  be  placed.  Also,  the  passive 
barrier  may  be  bypassed  due  to  fluctuations  in  groundwater  flow,  and  the 
effective  life  of  the  iron  filings  is  a  major  uncertainty  (NRC,  1999). 
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Another  strategy  that  has  gained  acceptance  over  the  years  as  a  low-cost  option 
is  natural  attenuation.  The  Environmental  Protection  Agency  defines  natural 
attenuation  in  the  following  manner: 


The  term  “monitored  natural  attenuation,”  as  used  in  this  Directive,  refers  to  the 
reliance  on  natural  attenuation  processes  (within  the  context  of  a  carefully 
controlled  and  monitored  clean-up  approach)  to  achieve  site-specific  remedial 
objectives  within  a  timeframe  that  is  reasonable  compared  to  that  offered  by 
other  more  active  methods.  The  “natural  attenuation  processes”  that  are  at  work 
in  such  a  remediation  approach  include  a  variety  of  physical,  chemical,  or 
biological  processes  that,  under  favorable  conditions,  act  without  human 
intervention  to  reduce  the  mass,  toxicity,  mobility,  volume,  or  concentration  of 
contaminants  in  soil  and  groundwater.  These  in-situ  processes  include 
biodegradation;  dispersion;  dilution;  sorption;  volatilization;  radioactive  decay; 
and  chemical  or  biological  stabilization,  transformation,  or  destruction  of 
contaminants  (USEPA,  1999). 


While  chlorinated  solvents  are  resistant  to  biodegradation,  they  have  been 
demonstrated  to  degrade  to  innocuous  compounds,  such  as  ethene,  under 
appropriate  conditions  (Maymo-Gatell  et.  al.,  1999).  Natural  attenuation  of 
chlorinated  compounds  in  an  aquifer  by  indigenous  microorganisms  is  typically 
very  cost  effective,  with  the  main  expenses  being  site  characterization  and 
monitoring.  However,  as  MNA  is  dependent  upon  natural  processes  taking  place 
in  the  subsurface  (as  compared  to  the  engineered  processes  discussed  above), 
the  contaminant-removal  mechanisms  are  the  least  understood.  Recently, 
natural  attenuation  has  come  under  criticism  for  being  used  in  situations  where 
its  effectiveness  has  not  been  adequately  demonstrated.  According  to  a  report 
from  the  National  Research  Council  (NRC),  “natural  attenuation  should  only  be 
accepted  as  a  formal  remedy  for  contamination  only  when  the  processes  are 
documented  to  be  working  and  are  sustainable”  (NRC,  2000).  The  EPA  states 
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natural  attenuation  “should  be  used  with  caution  commensurate  with  the 
uncertainties  associated  with  the  particular  applications,”  and  that  “the  hydrologic 
and  geochemical  conditions  favoring  significant  biodegradation  of  chlorinated 
solvents  sufficient  to  achieve  remediation  objectives  within  a  reasonable 
timeframe  are  anticipated  to  occur  only  in  limited  circumstances”  (USEPA,  1999). 
One  way  to  demonstrate  that  natural  attenuation  may  be  occurring  in  the 
subsurface  is  to  apply  a  model.  Several  computer  models  have  been  presented 
that  simulate  the  natural  attenuation  of  PCE  and  TCE  (Feng,  2000;  Clement  et 
al.,  1999),  but  few  have  been  validated  using  data  from  a  contaminated  site. 
Model  validation  could  result  in  increased  credibility  for  natural  attenuation  as  a 
remediation  strategy  for  chlorinated  solvent  contamination,  with  a  subsequent 
rise  in  implementation  and  reduction  in  overall  costs. 

The  purpose  of  this  research  is  to  validate  a  computer  model  by  comparing  its 
output  to  data  collected  from  the  field.  This  research  will  be  limited  to  PCE  and 
TCE  degradation  in  saturated  groundwater  systems.  A  review  of  current 
literature  will  focus  on  1 )  important  physico-chemical  and  biological  processes 
thought  to  take  place  in  the  subsurface,  2)  numerical  models  that  can  simulate 
those  processes,  3)  model  validation  strategies,  4)  previous  model  validation 
efforts,  and  5)  a  list  of  sites  that  could  support  future  validation  studies.  After  the 
literature  review,  a  model  will  be  selected  and  applied  to  a  chosen  site.  A 
simulation  will  be  run,  and  the  model’s  prediction  will  be  statistically  evaluated 
against  field  data.  It  is  hoped  that  this  method  of  model  validation  will  prove 
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useful  in  the  validation  of  other  contaminant  fate  and  transport  models.  It  is  also 
hoped  that  this  effort  will  generate  further  understanding  of  the  mechanisms  of 
natural  attenuation  of  chlorinated  solvents  in  aquifer  systems. 
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2.0  Literature  Review 


2.1  Overview 

The  purpose  of  this  literature  review  is  to  provide  background  information  on 
validation  methods  as  they  may  be  applied  to  modeling  natural  attenuation  of 
chlorinated  solvents  in  an  aquifer.  To  this  end,  this  review  covers  three  major 
areas.  First,  models  will  be  discussed,  including  the  definition  and  importance  of 
models  in  the  study  of  contaminated  aquifers,  the  natural  attenuation  processes 
thought  to  take  place  in  the  saturated  zone,  and  the  numerical  models  that 
simulate  natural  attenuation  processes.  Second,  the  concept  of  model  validation 
will  be  introduced,  and  will  include  a  discussion  of  model  verification,  calibration, 
and  comparative  analysis.  Third,  several  cases  involving  model  application  to 
chlorinated  ethene-contaminated  sites  will  be  reviewed  to  illustrate  common 
model  validation  practices. 

2.2  Models 

2.2.1  Definition,  Importance,  and  Uses 

A  model  may  be  defined  as  “a  representation  of  a  real  system  or  process” 
(Konikow  and  Bredehoeft,  1992).  Within  this  definition,  there  is  a  wide  array  of 
models  that  differ  greatly  in  purpose,  application,  structure,  and  complexity.  This 
research  will  focus  on  numerical  models  that  represent  the  fate  and  transport  of 
chlorinated  ethenes  in  an  aquifer. 
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Numerical  models  are  approximations  of  the  exact  mathematical  solution  of  the 
governing  equation(s)  that  describe  real  systems  and  processes.  These  models 
rely  on  fewer  simplifying  assumptions  than  analytical  models  that  solve  the 
governing  equation(s)  exactly,  and  are  capable  of  addressing  difficult  problems 
such  as  heterogeneous  conditions  and  complex  initial  and  boundary  conditions 
(Weaver,  et  al.,  1989).  The  Environmental  Protection  Agency  (EPA)  uses 
numerical  models  in  a  predictive  mode  to  make  regulatory  assessments  and 
environmental  decisions  (Weaver,  et  al.,  1989).  Numerical  models  have  also 
been  used  in  courtroom  litigation  to  establish  liability  (Bair,  1994).  Perhaps  more 
important  than  the  predictive  capability  of  a  model  is  its  explanatory  power. 
According  to  Murphy  and  Ginn  (2000),  who  specifically  comment  upon  modeling 
the  microbial  processes  that  occur  during  natural  attenuation,  “progress  in 
modeling  microbial  processes  in  porous  media  is  essential  to  improving  our 
understanding  of  how  physical,  chemical,  and  biological  processes  are  coupled 
in  groundwater  and  their  effect  on  groundwater-chemistry  evolution, 
bioremediation,  and  the  reactive  transport  of  contaminants  and  bacteria.”  The 
remainder  of  this  section  will  be  devoted  to  describing  the  processes  thought  to 
occur  in  chlorinated  ethene-contaminated  aquifers  and  the  numerical  models  that 
simulate  those  processes. 

2.2.2  Processes  Modeled 

Only  the  natural  attenuation  processes  thought  to  be  important  to  the  fate  and 
transport  of  chlorinated  ethenes  and  their  daughter  products  will  be  studied. 
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These  processes  can  be  divided  into  two  major  categories:  physiochemical 
processes  and  biological  processes. 

2.2.2. 1  Physicochemical  Processes 

The  physicochemical  processes  of  advection,  dispersion,  and  sorption  play  a 
crucial  role  in  the  fate  and  transport  of  contaminants.  Advection  is  the  transport 
of  mass  due  to  the  bulk  flow  of  groundwater,  and  is  “by  far  the  most  dominant 
mass  transport  processes”  (Dominico  and  Schwartz,  1998).  Dominico  and 
Schwartz  (1998)  provide  equations  describing  the  advective  transport  of 
contaminants. 

In  almost  all  cases,  mass  is  transported  beyond  the  region  delineated  solely  by 
advective  transport.  This  is  due  to  dispersion,  which  is  the  spreading  of  mass 
due  to  fluid  mixing.  Dispersion  is  caused  by  molecular  diffusion  (Brownian 
motion)  as  well  as  mechanical  mixing  within  a  heterogeneous  aquifer. 
Dispersion  is  often  modeled  as  diffusion,  as  the  outcome  of  each  mechanism  is 
similar  (Dominico  and  Schwartz,  1998).  Clark  (1996)  provides  a  detailed 
explanation  of  the  equations  used  to  describe  dispersion. 

Contaminants,  especially  non-polar  organic  compounds  such  as  PCE  and  TCE, 
typically  do  not  move  as  quickly  as  predicted  by  advection  and  dispersion.  This 
delayed  contaminant  movement,  or  retardation,  is  caused  by  several 
mechanisms  that  are  collectively  termed  sorption.  Sorption  is  defined  as  the 
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partitioning  of  contaminant  from  the  aqueous  phase  (i.e.  dissolved  in 
groundwater)  to  the  aquifer  solids.  In  this  research,  the  dominant  sorption 
process  involves  organic  contaminants  being  sorbed  to  the  organic  material 
found  on  aquifer  solids.  The  amount  of  sorption  that  takes  place  is  dependent 
upon  temperature,  concentration  of  the  contaminant,  and  the  characteristics  of 
both  the  contaminant  and  the  aquifer  solids. 

A  set  of  experiments  can  be  run  to  determine  the  relation  at  equilibrium  between 
aqueous  contaminant  concentration  (mass  of  contaminant  per  volume  water)  and 
sorbed  concentration  (mass  of  contaminant  per  mass  of  sorbent).  As  these 
experiments  are  performed  at  the  same  temperature,  the  relationship  is  referred 
to  as  a  sorption  isotherm. 

If  sorption  occurs  quickly  with  respect  to  groundwater  flow,  then  it  is  assumed 
that  the  rate  at  which  mass  is  being  sorbed  to  the  solids  is  equal  to  the  rate  at 
which  mass  is  being  desorbed.  This  is  called  the  local  equilibrium  assumption 
(LEA).  In  this  case,  the  partitioning  process  is  said  to  be  in  equilibrium  and  is  best 
described  as  an  equilibrium  sorption  isotherm.  Equilibrium  sorption  isotherms 
may  either  be  linear  or  non-linear.  According  to  Fetter  (1993),  a  linear  isotherm 
can  be  described  by  the  equation 

C*  =  Kd  C  (2.1) 

where 

C*  =  mass  of  contaminant  sorbed  per  dry  unit  weight  of  solid  (M/M) 
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C  =  concentration  of  contaminant  in  solution  (M/L3) 
Kd  =  distribution  coefficient  (L3/M) 

Figure  2.1  is  an  example  of  a  linear  isotherm. 


Figure  2.1  Linear  Sorption  Isotherm  (after  Fetter  (1993)) 

As  can  be  observed  from  Figure  2.1,  the  major  drawback  of  describing  sorption 
with  a  linear  isotherm  is  that  there  is  no  apparent  limit  to  the  amount  of 
contaminant  that  may  be  sorbed.  As  no  real  material  can  sorb  an  endless 
amount  of  mass,  linear  models  are  typically  used  only  at  relatively  low 
concentration  levels  or  within  a  small  range  of  contaminant  concentrations. 

When  describing  sorption  over  a  wide  range  of  contaminant  concentrations,  non¬ 
linear  isotherms  typically  perform  better  than  linear  isotherms.  Non-linear 
sorption  isotherms  describe  a  curvilinear  relationship  between  sorbed 
contaminant  concentration  and  dissolved  concentration.  The  Freundlich 
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isotherm  is  a  commonly  used  non-linear  sorption  isotherm,  and  is  described  by 
the  following  equation: 

C*  =  KCn  (2.2) 

where 

C*  =  mass  of  contaminant  sorbed  per  dry  unit  weight  of  solid  (M/M) 

C  =  concentration  of  contaminant  in  solution  (M/L3) 

K  =  constant 
N  =  constant  (-) 

Figure  2.2  illustrates  a  Freundlich  sorption  isotherm  with  K  =  28  and  N  =  0.62 
(Clark,  1996). 


Figure  2.2  Freundlich  Isotherm  (after  Fetter  (1993)) 

Note  that  the  Freundlich  isotherm  suffers  from  the  same  problem  as  the  linear 
isotherm.  That  is,  the  amount  of  contaminant  that  can  be  sorbed  by  aquifer 
solids  is  essentially  unlimited.  This  problem  can  be  addressed  through  the  use  of 
the  Langmuir  sorption  isotherm. 
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The  Langmuir  isotherm  is  a  non-linear  isotherm  that  is  based  on  the  assumption 
that  aquifer  solids  have  a  limited  number  of  sites  available  for  contaminant 
sorption.  The  Langmuir  isotherm  is  given  by 

C  =  1  C 

C*”a(3+(3  (2.3) 

where 

a  =  an  absorption  constant  related  to  the  binding  energy  (L3/M) 

P  =  the  maximum  amount  of  solute  that  can  be  absorbed  by  the  solid  (M/M) 

Figure  2.3  is  a  graph  of  a  Langmuir  isotherm  with  a  value  of  a  of  0.9  and  (3  of  0.9. 


Figure  2.3  Langmuir  Isotherm  (after  Fetter  (1993)) 

Not  all  groundwater  systems  can  be  described  with  equilibrium  isotherms.  In 
some  aquifers,  sorption  takes  place  at  approximately  the  same  rate  as  the 
velocity  of  groundwater.  In  these  systems,  local  equilibrium  cannot  be  assumed. 
Such  systems  are  more  appropriately  described  by  non-equilibrium,  or  rate- 
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limited,  sorption  models.  These  models  are  significantly  more  complicated  than 
equilibrium  isotherms,  and  are  beyond  the  scope  of  this  research.  Fetter  (1993) 
and  Feng  (2000)  contain  detailed  discussions  on  commonly  used  rate-limited 
models. 

2.2.2.2  Biological  Processes 

While  the  primary  physicochemical  processes  either  spread  the  contaminant 
(advection  and  dispersion)  or  temporarily  sequester  it  (sorption),  they  do  not 
reduce  the  overall  mass  of  contaminant  within  the  aquifer.  Conversely,  it  has 
been  demonstrated  that  microorganisms  are  capable  of  catalyzing  chemical 
reactions  that  result  in  the  degradation  of  chlorinated  ethenes  (Wilson  and 
Wilson,  1985).  This  research  will  cover  the  microbially-mediated  reductive  and 
oxidative  processes  known  to  degrade  chlorinated  ethenes  in  contaminated 
aquifers. 

Oxidation-reduction  (redox)  reactions  involve  the  transfer  of  an  electron  from  an 
electron-rich  chemical,  or  electron  donor,  to  an  electron-poor  chemical,  or 
electron  acceptor  (NRC,  2000).  Microorganisms  catalyze  redox  reactions  with 
enzymes  and  cofactors  in  order  to  generate  the  adenosine  triphosphate  (ATP) 
needed  to  sustain  metabolism  and  growth.  The  amount  of  ATP  generated 
depends  on  the  electron  donor  and  electron  acceptor  used  (NRC,  2000).  The 
most  common  oxidation-reduction  reactions  that  degrade  chlorinated  ethenes  are 
discussed  below. 


14 


2.2.2.2.1  Cometabolic  Reductive  Dehalogenation 

Some  reactions  do  not  contribute  to  microorganism  growth  or  metabolism.  Such 
reactions,  termed  cometabolic,  may  nevertheless  result  in  the  fortuitous 
degradation  of  contaminants.  The  cometabolic  reductive  dehalogenation  of 
chlorinated  ethenes  is  an  important  example  of  this  type  of  reaction.  In  this 
anaerobic  process,  a  non-specific  enzyme  acts  upon  a  chlorinated  ethene, 
replacing  its  chlorine  substituent  with  a  hydrogen  ion  and  two  electrons  (NRC, 
2000).  For  example,  PCE  (with  four  chlorine  substituents)  would  be  reduced  to 
TCE  (with  three  chlorine  substituents).  As  the  transformation  of  chlorinated 
ethenes  provide  no  energy  to  the  microorganisms,  metabolic  electron  acceptors 
such  as  nitrate  or  carbon  dioxide  are  required  for  ATP  production. 

Biodegradable  organic  materials,  such  as  natural  organic  matter  or  petroleum 
hydrocarbons,  act  as  the  necessary  electron  donors  (Wiedemeier,  1996). 

Cometabolic  reductive  dehalogenation  often  fails  to  completely  reduce  PCE  and 
TCE  to  ethene.  This  is  due  to  several  factors.  First,  the  chlorinated  solvents 
must  compete  with  the  metabolic  electron  acceptors  for  available  enzyme  sites. 
This  is  referred  to  a  competitive  inhibition.  As  the  solvents  only  gain  a  small 
share  of  the  available  electrons,  a  typical  aquifer  is  depleted  of  electron  donors 
well  before  the  contaminants  are  reduced  to  ethene.  Only  aquifers  co¬ 
contaminated  with  landfill  leachate  or  petroleum  hydrocarbons  see  significant 
reductive  dehalogenation.  Second,  the  solvents  or  their  daughter  products  may 
be  toxic  to  some  organisms,  thereby  negatively  affecting  some  microbe 
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populations  (Azadpour-Keeley  et  al.,  1999).  Third,  daughter  products  with  fewer 
chlorine  substituents  are  less  oxidized,  and  therefore  tend  to  be  reduced  less 
quickly,  if  at  all.  This  results  in  the  accumulation  of  cis-DCE  and  VC  (Major  et  al., 
1991 ;  Wilson  et  al.,  1995),  which  is  especially  problematic  in  that  VC  is  more 
toxic  than  PCE  orTCE  (Masters,  1997).  For  these  reasons,  cometabolic 
reductive  dehalogenation  is  “considered  ubiquitous  in  anaerobic  systems  but 
generally  incapable  of  mediating  complete  reduction  to  non-toxic  products  like 
ethene”  (Bradley,  2000). 

2.2.2.2.2  Reductive  Dehalogenation  Through  Halorespiration 

It  has  been  recently  discovered  that  not  all  reductive  dechlorination  processes 
are  cometabolic.  A  group  of  microorganisms,  called  halorespirers,  have  been 
found  to  generate  ATP  by  using  chlorinated  ethenes  as  sole  terminal  electron 
acceptors  (Holliger  et  al.,  1993;  Maymo-Gatell  et  al.,  1997).  As  these 
microorganisms  derive  energy  from  the  process,  under  the  appropriate 
conditions  they  are  capable  of  higher  rates  of  dechlorination  than  cometabolizers 
(Bradley,  2000).  Indeed,  according  to  Wiedemeier  et  al.  (1999),  halorespiration 
“probably  accounts  for  the  majority  of  chlorinated  solvent  biodegradation  at  many 
of  the  sites  where  biodegradation  is  significantly  attenuating  the  [chlorinated 
solvent]  plume.”  However,  the  electron  donors  that  are  used  by  halorespirers  are 
limited  to  hydrogen  and  possibly  acetate  and  formate.  As  sulfate  reducers, 
methanogens,  and  homoacetogens  also  compete  for  these  electron  donors,  the 
ability  of  halorespirers  to  reduce  VC  to  ethene  may  be  significantly  diminished 
(McCarty,  1996;  Smatlak  et  al.,  1996). 
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2. 2. 2. 2.3 


Cometabolic  Oxidation 


As  more  chlorine  substituents  are  removed  and  the  compound  becomes  more 
reduced,  the  chlorinated  ethene  may  be  more  easily  oxidized.  A  number  of 
microorganisms  have  been  identified  that  can  oxidize  TCE  and  its  daughters  (but 
not  PCE)  to  CO2  cometabolically  (McCarty  and  Semprini,  1994).  This  may  be  an 
important  process  when  oxygen  and  the  necessary  carbon  substrate  (e.g. 
methane,  ethylene,  phenol,  toluene)  coexist,  such  as  at  the  fringe  of  a 
contaminant  plume  (Dolan  and  McCarty,  1995;  Anderson  and  McCarty,  1997). 

2.2.2.2.4  Aerobic  Oxidation 

In  some  cases,  DCE  and  VC  may  be  oxidized  directly  in  a  process  that  is 
beneficial  to  the  microorganisms.  In  aerobic  oxidation,  VC  can  be  used  as  a  sole 
carbon  source  for  growth  and  metabolism,  while  DCE  has  been  shown  as  a 
carbon  substrate  for  metabolism  only.  As  with  cometabolic  oxidation,  aerobic 
oxidation  may  be  important  at  the  fringe  of  a  plume  (Bradley  and  Chapelle, 
1998b). 

2.2.2.2.5  Anaerobic  Oxidation 

Anaerobic  oxidation  is  another  possible  pathway  for  chlorinated  ethene 
degradation.  Microorganisms  have  been  shown  to  oxidize  VC  to  CO2  under  Fe 
(III)— reducing  conditions.  DCE  has  been  shown  to  oxidize  directly  to  C02  under 
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Mn  (IV)-reducing  conditions  (Bradley  et  al.,  1998).  As  these  processes  can  take 
place  in  an  anaerobic  environment  in  conjunction  with  reductive  dehalogenation, 
this  may  prove  to  be  a  significant  pathway  for  chlorinated  ethene  degradation. 


2.2.2.2.6  Mathematical  Description  of  Contaminant  Biodegradation 

If  there  is  no  shortage  of  necessary  substrates  (i.e.  electron  donors  and 
acceptors)  then  contaminant  degradation  can  be  expressed  as  a  first  order 
decay.  That  is,  the  rate  at  which  the  contaminant  disappears  is  dependent  only 
upon  the  contaminant  concentration.  This  is  modeled  by  the  equation 


dt 


(2.4) 


where 

C  is  the  concentration  of  dissolved  contaminant  (M/L3) 
k  is  the  contaminant  decay  first  order  rate  constant  (T'1) 
Figure  2.4  is  an  example  of  a  first  order  process. 


t  (sec) 


Figure  2.4  First  Order  Biodegradation 
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Reductive  dehalogenation  is  typically  modeled  as  a  first  order  process,  with  a 
separate  decay  constant  for  each  reductive  step  (Clement  et  al.,  2000).  The 
direct  oxidation  of  chlorinated  hydrocarbons  is  also  modeled  as  first  order 
(Bradley  and  Chapelle,  1998a). 


In  situations  where  one  or  more  substrates  limit  the  rate  of  biodegradation,  it  is 
necessary  to  use  Monod  kinetics  to  describe  the  rate  of  the  reaction.  Monod 
kinetics  describe  the  growth  of  microorganisms  on  a  limiting  substrate  (Suarez 
and  Rifai,  1999)  by  the  hyperbolic  saturation  function 


M- 


M' 


max 


f  S  1 

ls  +  KsJ 


(2.5) 


where 

p.  =  microorganism  growth  rate  (T'1) 

Pmax  =  maximum  growth  rate  of  microorganisms  (T'1) 
Ks  =  half  saturation  constant  (M/L3) 

S  =  concentration  of  limiting  substrate  (M/L3) 


The  half  saturation  constant  is  the  concentration  of  limiting  substrate  at  which  the 
microorganisms  grow  at  half  the  maximum  growth  rate  (Wiedemeier  et  al.,  1999). 
Dual  -  Monod  kinetics  have  been  used  to  describe  aerobic  cometabolism 
(Bouwer  and  McCarty,  1985),  as  the  concentration  of  both  electron  acceptors 
and  contaminants  limit  the  rate  of  biodegradation.  This  model  is  expressed  as 
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(2.6) 


— — x-kl 

f  C  1 

f  C*  1 

dt 

^Kg  +  Cy 

,KSA  +Ca; 

where 

X  =  concentration  of  microorganisms  active  for  cometabolism  (M/L3) 
k  =  maximum  utilization  rate  of  cometabolism  (M/M  -T'1) 

C  =  concentration  of  target  contaminant  (M/L3) 

Ks  =  half  saturation  constant  of  the  target  contaminant  (M/L3) 

Ca  =  concentration  of  electron  acceptor  (M/L3) 

Ksa  =  half  saturation  constant  of  electron  acceptor  (M/L3) 


The  above  model  does  not  take  into  consideration  the  competitive  inhibition 
between  the  electron  donor  and  the  contaminant  for  available  enzyme  sites.  The 
following  modification  of  (2.6)  accounts  for  this  competition  (Semprini  and 
McCarty,  1992): 


f 


\ 


dC 

dt 


=  -X  ■  k 


Ka.  CD 
Ko  +  C  h — - 


K 


sD  ) 


kSa+ca 


(2.7) 


where 

Cd  =  electron  donor  concentration  (M/L3) 

KSd  =  half  saturation  constant  of  electron  donor  (M/L3) 

Note  that  X,  the  concentration  of  microorganisms  found  in  equations  (2.6)  and 
(2.7),  also  depends  upon  electron  donor  and  acceptor  concentrations.  Monod 
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kinetics  can  be  used  to  describe  this  relationship  as  well,  as  in  the  equation 
proposed  by  Semprini  and  McCarty  (1991),  where  dual-Monod  kinetics  describes 
growth  and  Monod  kinetics  describes  decay: 


dX 

dt 


=  -XkY 


Y 


Ksd  +CDy 


Kc  +  C, 


bX 


ksa+Ca 


where 

k  =  maximum  utilization  rate  of  cometabolism  (M/M  -  T'1) 

Y  =  yield  coefficient  of  biomass  produced  per  substrate  used  (M/M) 
b  =  microbial  decay  rate  constant  (T'1) 


(2.8) 


It  is  important  to  note  that  the  above  equations  are  general,  and  that  the  set  of 
parameter  values  used  will  depend  upon  the  specific  oxidation-reduction  reaction 
being  modeled.  For  example,  the  parameters  used  to  model  the  reductive 
dehalogenation  of  PCE  to  TCE  with  sulfate  as  an  electron  acceptor  will  differ 
from  the  parameters  used  to  model  the  reduction  of  TCE  to  DCE  under 
methanogenic  conditions.  Therefore,  a  model  that  simulates  all  relevant 
processes  should  be  able  to  track  the  concentration  of  multiple  electron  donors, 
acceptors,  and  contaminants,  then  apply  the  appropriate  parameters  as  required. 
Typically,  simplifying  assumptions  are  built  into  models,  and  multiple  donors  and 
acceptors  are  not  tracked.  Of  the  models  discussed  in  the  next  section,  only  Bio- 
Redox  accounts  for  reactions  among  multiple  electron  donors  and  acceptors. 
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2.2.3  Models  Relevant  to  the  Simulation  of  Chlorinated  Ethene  Natural 
Attenuation 

There  are  many  models  available  that  may  be  used  to  simulate  chlorinated 
ethene  fate  and  transport.  To  reduce  the  number  of  models  to  be  considered  in 
detail,  an  initial  screening  was  accomplished  using  several  criteria.  Due  to  the 
time  and  money  available  for  this  research,  only  completed  models  that  were 
readily  available  at  little  or  no  cost  were  considered.  As  discussed  in  Chapter  1 , 
this  research  focuses  on  the  natural  attenuation  of  aqueous  phase  chlorinated 
solvents  in  an  aquifer.  Accordingly,  models  that  simulated  transport  in  the 
unsaturated  zone,  a  combination  of  the  saturated  and  unsaturated  zone,  or 
simulated  dual-phase  flow  were  rejected.  Finally,  only  those  models  that  were 
able  to  simulate  the  reactive  step-wise  degradation  of  chlorinated  ethenes  were 
investigated.  After  this  initial  screening,  four  models  were  examined  in  some 
detail:  Biochlor,  Bio-Redox,  RT3D,  and  BR3D. 

2.2.3.1  Biochlor 

Biochlor  is  described  as  a  "natural  attenuation  decision  support  system"  (Aziz  et 
al.,  2000)  for  sites  with  dissolved  chlorinated  solvents.  This  computer  code 
provides  concentration  data  along  the  plume  centerline  by  describing  one¬ 
dimensional  (1-D)  advection,  three-dimensional  (3-D)  dispersion,  linear  sorption, 
and  biotransformation  due  to  reductive  dehalogenation  (Aziz  et  al.,  2000).  The 
biotransformation  process  is  modeled  as  a  sequential  first  order  decay,  and  two 
separate  reaction  zones  can  be  simulated. 
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Biochlor  is  a  screening  model  that  assumes  simple  groundwater  flow,  uniform 
hydrogeologic  and  environmental  conditions,  and  a  constant  vertical  plane 
source  (Aziz  et  al.,  2000).  The  purpose  of  this  model  is  to  facilitate  remediation 
planning,  and  it  is  not  designed  to  accurately  predict  contaminant  concentrations 
in  a  complex,  real  aquifer. 

2.2.3.2  Bio-Redox 

Bio-Redox  is  a  comprehensive  3-D  model  used  to  simulate  multiple  oxidation- 
reduction  reactions  (Feng,  2000).  The  model  describes  3-D  advection  and 
dispersion,  linear  and  non-linear  equilibrium  sorption,  but  not  rate-limited 
sorption.  Bio-Redox  can  simulate  contaminant  biotransformation  due  to 
reductive  dehalogenation,  direct  oxidation,  or  methanotrophic  cometabolism  by 
using  either  first  order,  Monod,  Dual-Monod  kinetics  (with  and  without 
competitive  inhibition).  This  computer  code  can  simulate  reactions  among 
multiple  electron  donors  and  acceptors,  and  can  track  the  accumulation  of 
chloride  ions  released  by  reductive  dehalogenation.  As  the  application  of  Bio- 
Redox  to  a  site  does  not  require  any  modification  of  the  FORTRAN  source  code, 
the  program  is  relatively  simple  to  use. 

2.2.3.3  RT3D 

Reactive  Multi-species  Transport  in  3-Dimensions,  otherwise  known  as  RT3D, 
simulates  multiple  transport  equations  coupled  with  multiple  biochemical  kinetics 
(Clement,  2000).  RT3D  is  a  modular  program,  composed  of  subprograms  that 
allow  the  modeler  to  describe  3-D  advection  and  dispersion;  sorption  through 
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linear  equilibrium,  non-linear  equilibrium,  or  rate  limited  reactions;  and 
biodegradation  through  first  order,  Monod,  or  Dual-Monod  kinetics.  The  modular 
design  of  RT3D  makes  it  highly  adaptable  to  a  specific  site;  seven 
preprogrammed  modules  are  available,  and  modelers  can  create  their  own 
module  to  run  within  the  program.  While  RT3D’s  design  allows  for  great  flexibility 
and  creativity  in  applying  the  model  to  a  site,  the  level  of  effort  required  is 
considerable.  The  appropriate  reaction  module  must  be  chosen  (Clement,  1997) 
and  correctly  placed  within  the  FORTRAN  source  code,  a  task  that  is  difficult  due 
to  the  code’s  complexity  (Feng,  2000).  The  unintentional  modification  of  source 
code  is  also  of  concern. 

2.2.3.4  BR3D 

BR3D  is  a  model  developed  at  the  Air  Force  Institute  of  Technology  (AFIT)  to 
simulate  the  degradation  of  chlorinated  ethenes  through  the  cometabolic 
processes  of  reductive  dechlorination  and  aerobic  oxidation.  The  capabilities  of 
BR3D  are  similar  to  RT3D  in  that  it  can  simulate  3-D  advection  and  dispersion; 
sorption  through  linear  equilibrium,  non-linear  equilibrium,  or  rate  limited 
reactions;  and  biodegradation  through  first  order,  Monod,  or  Dual  Monod  kinetics. 
It  cannot  combine  multiple  electron  acceptors  and  donors  in  oxidation-reduction 
processes  (Feng,  2000).  BR3D  is  relatively  easy  to  use,  as  it  does  not  require 
modification  of  the  FORTRAN  source  code. 
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2.3  Model  Validation 


There  is  no  single  definition  of  the  term  “validation”  within  the  scientific 
community  (Leijnse  and  Hassanizadeh,  1994).  The  International  Atomic  Energy 
Agency,  an  entity  that  supervises  the  development  of  many  long-term  predictive 
models,  defines  validation  in  the  following  manner  (IAEA  1982): 

A  conceptual  model  and  the  computer  code  derived  from  it  are  validated  when  it 
is  confirmed  that  the  conceptual  model  and  the  computer  code  provide  a  good 
representation  of  the  actual  processes  occurring  in  the  real  system. 

Schlesinger  (1979)  provides  a  similar  definition,  stating  that  validation  is 

Substantiation  that  a  computerized  model  within  its  domain  of  applicability 
possesses  a  satisfactory  range  of  accuracy  consistent  with  the  intended 
application  of  the  model. 

A  considerable  number  of  arguments  about  model  validation  are  centered  on  the 
scope  of  the  validation  effort.  Sargent  (1982)  states  that  validation  consists  of 
three  components:  validation  of  the  logic  and  consistency  of  the  model 
(conceptual  validation);  determination  of  the  model’s  ability  to  answer  the 
question  at  hand  (operational  validation);  and  comparison  of  the  simulation  to  the 
observed  system  (quantitative  validation).  Others  argue  that  validation  only 
consists  of  the  comparison  of  model  output  to  independent  observations 
(McCombie  and  McKinley,  1998;  ASTM,  1996;  USEPA,  1989).  For  the  purposes 
of  this  work,  the  overarching  concept  of  building  confidence  in  a  model’s 
predictive  capability  will  be  referred  to  as  model  validation,  while  the  specific  task 
of  comparing  model  output  to  independent  observations  will  be  termed 
comparative  analysis. 
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In  a  strict  sense,  no  model  can  ever  be  considered  true  (i.e.  valid)  under  all 
circumstances  (Tsang,  1991).  In  accordance  with  the  scientific  method,  a 
hypothesis  can  be  proved  ‘false’  or  ‘not  false,’  but  may  never  be  considered 
‘true,’  ‘valid,’  or  ‘substantiated’  (Popper,  1959;  Konikow,  1992;  Bredehoeft  and 
Konikow,  1 993).  However,  Niederer  argues  that  “it  does  not  make  sense  to 
demand  strict  proof  that  a  model  is  correct;  it  makes  a  lot  of  sense,  however,  to 
promote  consensus  by  providing  ample  positive  evidence  for  the  correctness  of 
the  model.  In  this  sense,  validation  is  primarily  a  means  to  achieve  consensus” 
(Niederer,  1990).  In  order  to  build  consensus  that  a  model  accurately  describes 
important  processes,  the  structure  of  the  model  itself  should  be  verified,  and  its 
application  to  a  particular  problem  should  be  validated  through  model  calibration 
and  comparative  analysis  of  model  output  to  observed  conditions. 

2.3.1  Model  Verification 

A  model’s  structure  is  verified  through  an  extensive  process  of  documentation 
and  testing  of  the  computer  code.  Documentation  of  a  computer  code  begins 
with  the  development  of  the  algorithms  and  procedures,  and  continues 
throughout  the  modeling  process  (ASTM,  1996).  This  allows  others  to  check  the 
accuracy  of  the  code  as  well  as  understand  the  assumptions  upon  which  the 
model  is  based.  Computer  code  verification  is  “the  process  of  demonstrating  the 
consistency,  completeness,  correctness,  and  accuracy  of  a  ground-water 
modeling  code  with  respect  to  its  design  criteria  by  evaluating  the  functionality 
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and  operational  characteristics  of  the  code  and  testing  embedded  algorithms  and 
internal  data  transfers  through  execution  of  problems  for  which  independent 
benchmarks  are  available”  (ASTM,  1996).  Analytical  results  or  results  from  a 
known  data  set  are  common  independent  benchmarks  used  to  ensure  that  the 
computer  code  is  performing  as  expected  (Clement  et  al.,  2000).  A  computer 
code  is  verified  when  it  is  determined  to  be  “mathematically  correct  in  the 
formulation  and  solution”  (Tsang,  1991). 

2.3.2  Model  Calibration 

Model  verification  is  no  guarantee  that  the  assumptions  and  processes  simulated 
in  a  model  are  applicable  to  a  given  real  system.  To  show  that  the  model  can 
represent  such  systems  requires  model  calibration.  Calibration  is  the  process  of 
adjusting  model  parameters,  initial  and  boundaries  conditions,  and  stresses  so 
that  the  model  approximates  field-measured  values  (Levy,  1993).  Calibration  is 
an  important  step  in  model  validation,  as  it  can  uncover  flaws  in  the  design  and 
implementation  of  the  model.  According  to  the  NRC,  “if  the  model  cannot 
capture  the  observed  trends  no  matter  how  well  it  is  calibrated,  then  its 
conceptual  basis  surely  is  wrong”  (NRC,  2000).  Calibration  is  particularly  useful 
when  necessary  parameter  values  are  unavailable  or  difficult  to  obtain.  The 
standard  calibration  procedure  seeks  to  minimize  the  differences  between 
observed  and  predicted  behavior,  as  measured  by  some  goodness-of-fit  statistic 
(Armstrong  et  al.,  1996).  Parameter  values  that  optimize  the  goodness-of-fit 
statistic  can  be  determined  through  the  use  of  nonlinear  regression  analysis  (Hill, 
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1998),  or  through  a  trial-and-error,  forward  process  (Levy,  1993;  Resele  and  Job, 
1990).  Mean  squared  error  (MSE)  is  a  commonly  used  goodness-of-fit  statistic, 
and  will  be  discussed  in  later  sections.  Hill  (1998)  uses  nonlinear  regression  in 
the  computer  program  UCODE  to  optimize  the  parameter  values  (P’P)  by 
minimizing  the  following  statistic: 

ND  NPR 

S(b)  =  2>i[Qoi  -  Qsi]2  +  2>p[Pp  -  P'p  ]2  (2.9) 

i=1  p=1 

where 

b  =  a  vector  with  values  for  each  of  the  NP  parameters  being  estimated 
ND  =  number  of  observations 

NPR  =  number  of  independently  estimated  parameter  values 
NP  =  number  of  estimated  parameters  (NP>NPR) 

Qoi  =  ith  observation 

Qsi  (b)  =  simulated  value  corresponding  to  the  i,h  observation 
pp  =  pth  independently  estimated  parameter  value 
pp  ’  =  p,h  fitted  parameter  value 

cos  =  weight  for  the  ith  observation 

cop  =  weight  for  the  pth  independently  estimated  parameter 

Weighting  performs  two  functions  (Hill,  1998).  First,  the  weights  ensure  that  the 
weighted  residuals  have  the  same  units  so  that  they  may  be  squared  and 
summed  as  shown  in  equation  (2.9).  Second,  the  weights  determine  the  amount 
of  influence  an  observation  or  independently  estimated  parameter  value  has.  As 
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the  weights  are  based  on  the  uncertainty  of  the  measurement,  more  accurate 
estimates  and  observations  will  have  greater  influence  on  the  overall  statistic. 

The  benefit  of  this  statistic  is  that  it  penalizes  those  models  whose  fitted 
parameter  values  differ  significantly  from  independently  estimated  values. 
Parameter  values  should  be  chosen  carefully  if  such  a  penalty  is  not  incurred,  as 
it  is  possible  to  obtain  an  excellent  fit  between  predicted  and  observed  results 
using  unrealistic  parameter  values  (Bobba,  1993).  To  prevent  this,  parameter 
values  used  in  a  model  are  typically  bounded  to  values  that  are  deemed  realistic 
based  on  sound  judgment  (Bobba,  1993)  or  by  the  literature  (Clement  et  al., 
2000). 

Some  models  require  no  calibration,  as  they  are  used  in  a  purely  predictive 
manner  (Eggleston,  2000;  Bond,  1998;  Armstrong,  1996).  That  is,  the  input 
parameter  values  are  derived  independently  from  field  and  lab  experiments  or 
from  the  relevant  literature.  These  models  tend  to  be  relatively  simple  models, 
such  as  pesticide  leaching  models  or  simplified  groundwater  flow  models,  and 
can  be  difficult  to  use  when  parameters  values  cannot  be  easily  obtained 
(Clement,  2000).  According  to  the  NRC  (2000),  any  comprehensive  reactive 
transport  model  must  be  calibrated. 

Not  all  parameter  values  used  in  a  comprehensive  model  have  to  be  determined 
exclusively  by  calibration.  Through  the  use  of  geostatistical  methods  such  as 
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kriging,  parameter  values  such  as  hydraulic  conductivity  may  be  estimated.  This 
method  can  provide  a  more  independent  representation  of  the  real  system. 

2.3.3  Comparative  Analysis 

Once  a  model  has  been  calibrated,  simulations  can  be  run  and  compared  to 
additional  laboratory  or  field  data  outside  the  calibration  set.  According  to  the 
EPA  (1989),  “field  evaluations  of  models  often  lead  to  a  better  understanding  of 
the  processes  taking  place  and  point  to  additional  research  needs.”  This 
research  will  compare  model  results  to  field  observations.  Comparing  model 
results  to  field  observations  has  been  referred  to  as  groundtruth  model  testing 
(Ababou  et  al.,  1992)  or  model  validation  (McCombie  and  McKinley,  1998; 

ASTM,  1996;  USEPA,  1989).  For  the  purposes  of  this  research,  such  an 
evaluation  will  be  called  a  comparative  analysis.  Comparative  analysis  can  be 
divided  into  two  categories:  qualitative  comparative  analysis  and  quantitative 
comparative  analysis. 

2.3.3.1  Qualitative  Comparative  Analysis 

Qualitative  comparative  analysis  is  a  common  means  to  build  consensus  that  a 
model  is  simulating  the  processes  taking  place  in  the  field.  In  the  realm  of 
groundwater  contaminant  fate  and  transport  modeling,  model  output  is  typically 
compared  to  field  observations  through  the  use  of  concentration  plots.  Examples 
of  such  comparisons  include  contaminant  breakthrough  curves  and  two- 
dimensional  concentration  contour  plots.  While  this  method  allows  for  a  concise 
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and  informative  study  of  the  model’s  performance,  it  has  its  drawbacks.  Most 
notably,  qualitative  methods  are  subjective  in  nature,  which  makes  it  difficult  to 
compare  results  between  two  or  more  models.  To  overcome  this  drawback,  it  is 
recommended  that  qualitative  methods  be  combined  with  quantitative  methods  to 
produce  the  most  comprehensive  description  of  model  performance  (Willmott, 
1984;  Legates  and  McCabe,  1999). 

2.3.3.2  Quantitative  Comparative  Analysis 

Just  as  there  is  no  one  definition  of  validation,  there  is  also  no  one  best  statistic 
to  determine  the  goodness-of-fit  between  observations  and  model  output 
(Weglarczyk,  1998;  Imam  et  al.,  1998).  Accordingly,  various  combinations  of 
statistical  measures  are  used  to  evaluate  models  (Legates  and  McCabe,  1999; 
Imam  et  al.,  1998).  One  of  the  main  purposes  of  this  research  is  to  choose  a 
combination  of  statistics  that  can  be  used  to  evaluate  the  goodness-of-fit  of  a 
model  simulation.  This  section  will  discuss  the  residual-based  and  association- 
based  statistics  that  may  be  suitable  to  measure  model  goodness-of-fit. 

One  method  of  quantifying  model  performance  is  to  calculate  the  summation  of 
all  the  differences,  or  residuals,  between  the  simulated  values  (Qs)  and  their 
corresponding  observed  values  (0o).  Common  residual-based  statistics  include 
bias  (B),  sum  of  squares  error  (SSE),  and  mean  squared  error  (MSE). 

Bias  is  a  measure  of  the  systematic  error  associated  with  a  model  (Devore, 

1995).  The  amount  that  a  model  consistently  overestimates  or  underestimates 
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the  values  of  interest  is  reflected  in  the  calculation  of  bias.  Figure  2.5  illustrates 
the  difference  between  biased  and  unbiased  models.  Unconditional  bias  (B),  is 
defined  as 


B  =  m  s  -  m  o 


(2.10) 


where 

ms  =  average  simulated  value 
m0  =  average  observed  value 


-e —  Simulated 
Values 

■  Observed 
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Figure  2.5  Example  of  a  Model  a)  with  Bias,  and  b)  without  Bias 

If  a  nondimensional  measure  of  bias  is  preferred,  B'2  can  be  used  (Weglarczyk, 
1998).  B/2  is  defined  as 


(2.11) 
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where 


So  =  standard  deviation  of  the  observed  values 

Conditional  bias  (C2),  another  measure  of  model  performance,  is  the  measure  of 
covariance  between  model  residuals  (Qs  -  Qo,  or  AQ)  and  simulated  values  (Qs), 
and  is  defined  as 

C2  =  cov(Qs  -Q0,Qs)  =  s|  -Q0QS  -m0ms  (2.12) 

where 

Ss  =  standard  deviation  of  the  simulated  values 

Note  that  the  overline  indicates  the  average  of  the  quantities  indicated. 
Covariance,  an  association-based  statistic,  will  be  discussed  in  detail  later  in  this 
section.  A  completely  unbiased  model  will  have  values  of  B,  B'2,  and  C2  equal  to 
zero.  The  relationship  between  the  B/2,  and  C2,  and  other  statistics  will  also  be 
discussed  in  later  paragraphs. 

A  model  may  be  perfectly  unbiased  but  still  be  unacceptable  in  terms  of  matching 
simulated  values  to  observed  values.  Large  overestimates  may  be  balanced  by 
large  underestimates,  as  shown  in  Figure  2.6.  Statistics  that  sum  the  square  of 
errors  or  the  absolute  errors  are  required  to  better  assess  overall  model 
performance.  Calculations  of  this  nature  include  sum  of  squares  error  (SSE), 
mean  absolute  error  (MAE),  and  mean  squared  error  (MSE). 
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Time  (sec) 

Figure  2.6  Unbiased  Model  that  Fails  to  Match  Observed  Data 

Sum  of  squares  error  is  defined  as 

SSE=(Q0-Qs)2  (2.13) 

Mean  absolute  error  is  another  familiar  residual-based  statistic,  and  is  defined  as 

MAE  =  |Q0  -Qs|  (2.14) 

Even  more  commonly  used  is  the  statistic  mean  squared  error,  which  is  shown 
as 

MSE  =  (Q0  -  Qs  )2  (2.15) 

Root  mean  square  error  (RMSE),  the  square  root  of  MSE,  is  also  used  frequently 
as  it  has  the  same  units  as  the  residuals.  RMSE  has  been  used  extensively  as 
an  objective  function  in  model  calibration  (Sorooshian  et  al.,  1983),  while  MSE  is 
commonly  used  in  meteorological  forecast  validation  (Wilks,  1995).  Fox  (1981) 
and  Willmott  (1982)  consider  RMSE  to  be  among  the  best  residual-based 
performance  measures.  However,  as  both  MSE  and  RMSE  are  calculated  with 
squared  differences,  they  are  overly  sensitive  to  extreme  values  (Legates  and 
McCabe,  1999).  Additionally,  MSE  and  RMSE  are  cumbersome  to  use  when 
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comparing  the  performance  of  different  models,  as  their  values  are  dependant 
upon  the  units  of  the  observations  and  simulations.  As  a  result  of  this 
dimensionality,  the  MSE  of  a  model  that  is  simulating  data  measured  in 
milligrams  per  liter  (mg/L)  would  differ  significantly  from  a  model  that  simulates 
data  measured  in  moles  per  liter  (mol/L).  Non-dimensional  statistics  (discussed 
later  in  this  section)  should  be  used  along  with  MSE  and  RMSE  in  order  to 
quantify  model  goodness-of-fit  (Legates  and  McCabe,  1999). 

Residual-based  statistics  such  as  MSE  are  often  used  in  conjunction  with 
association-based  statistics  to  provide  the  modeler  with  a  comprehensive 
measure  of  model  quality.  To  best  understand  association-based  statistics, 
consider  the  scatter  plot  in  Figure  2.7.  Each  point  represents  a  measured  and 
simulated  value  at  a  certain  location  or  time.  The  x-coordinate  of  each  point  is 
determined  by  the  observed  value  (Q0),  while  the  simulated  value  (Qs) 
determines  the  y-coordinate.  Association-based  statistics  quantify  the 
correlation,  or  strength  of  the  linear  relationship,  between  the  predicted  and 
observed  values.  It  should  be  pointed  out  that  this  is  not  the  same  as  regression, 
which  implies  a  causal  relationship  between  independent  and  dependent 
variables.  Covariance  (cov),  sample  correlation  coefficient  (r),  coefficient  of 
determination  (R2),  and  coefficient  of  efficiency  (E),  are  commonly  used 
association-based  statistics. 
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Figure  2.7  Plot  of  Simulated  and  Observed  Values 

Covariance  measures  the  strength  of  the  relationship  between  the  observed  and 
simulated  values,  and  is  defined  as 

cov(Q0,Qs)  =  Q0QS  - m0ms  (2.16) 

where 

Qo  =  observed  values 
Qs  =  simulated  values 
ms  =  average  simulated  value 
m0  =  average  observed  value 

As  with  MSE,  covariance  is  unwieldy  due  to  its  dimensionality.  To  overcome  this, 
the  nondimensional  measure  of  covariation,  r,  can  be  used,  r,  known  as  the 
Pearson's  product-moment  correlation  coefficient  (Legates  and  McCabe,  1999), 
is  defined  by  Weglarczyk  (1998)  as 
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(2.17) 


_  cov(Q0,Qs)  _  Q0QS  m0ms 

S0Ss  S0Ss 

where 

So  =  standard  deviation  of  the  observed  values 
Ss  =  standard  deviation  of  the  simulated  values 

The  Pearson's  correlation  coefficient  is  bounded  from  -1  (the  largest  negative 
correlation)  to  1  (the  largest  positive  correlation).  A  more  common  measure  of 
correlation  is  R2,  the  coefficient  of  determination.  R2  is  calculated  as  the  square 
of  the  sample  correlation  index  (R2  =  r2)  and  describes  the  total  variance  in  the 
observed  data  that  can  be  explained  by  a  linear  model  of  correlation.  A  value  of 
R2  =1  is  achieved  when  all  points  lie  along  a  straight  line  (Devore,  1995).  It  is 
important  to  note  that  R2  “estimates  the  concentration  of  [(Qo,  Qs)]  points  along 
an  arbitrary  line  on  the  [(Q0,Qs)]  plane,  not  along  the  1 :1  line  which  is  of  the  only 
interest  to  the  modeller”  (Weglarczyk,  1998).  This  point  is  illustrated  in  Figure 
2.8,  where  two  modeling  results  are  displayed.  The  value  of  R2  is  the  same  for 
both  models,  but  Model  2  clearly  overestimates  the  observed  results,  as  denoted 
by  the  1 :1  line.  Due  to  the  insensitivity  of  r  and  R2  to  bias,  it  has  been  argued 
that  neither  should  be  used  as  a  measure  of  model  performance  (Legates  and 
McCabe,  1999;  Fox,  1981;  Willmott,  1981;  1984). 
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Figure  2.8  Insensitivity  of  Sample  Correlation  Coefficient  to  Bias 

One  association-based  statistic  that  does  account  for  model  bias  is  the  Nash- 
Sutcliffe  coefficient  of  efficiency,  E.  The  coefficient  of  efficiency  is  an  estimate  of 
the  concentration  of  (Qo,Qs)  points  along  the  1 :1  line,  or  line-of-perfect-fit.  The 
coefficient  of  efficiency  is  also  related  to  the  coefficient  of  determination  by  the 
relationship 

E  =  R2  -  C2  -  B'2  (2.18) 

The  coefficient  of  efficiency  can  also  be  calculated  as  a  dimensionless 
transformation  of  MSE,  and  is  defined  by 


E  =  1 


MSE 


(2.19) 


The  Nash-Sutcliffe  coefficient  has  several  desirable  characteristics.  First,  this 
coefficient  is  nondimensional,  making  it  easier  to  compare  the  performance  of 
different  models.  Second,  E  is  a  measure  of  the  departure  of  the  (Qo,Qs)  values 
from  the  line-of-perfect-fit,  which  is  highly  relevant  to  modelers.  Finally,  the 
coefficient  of  efficiency  increases  as  model  goodness-of-fit  increases.  A 
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maximum  value  of  E  =  1  suggests  perfect  model  performance,  while  a  negative 
value  of  E  indicates  that  the  model  “introduces  more  ambiguity  than  that 
introduced  by  simply  using  the  mean  value  of  the  observation  as  an  estimator” 
(Imam  et  al.,  1999).  While  E  is  a  useful  goodness-of-fit  statistic,  it  is  more  difficult 
to  interpret  than  R2.  According  to  Legates  and  McCabe  (1999),  an  E  of  0.70 
means  that  the  mean  square  error  accounts  for  30%  of  the  variance  in  the 
observed  data.  Frankenberger  et  al.  (1999)  describe  E  in  terms  of  mean 
observation  values,  such  that  an  E  of  0.70  indicates  that  the  model  performs  70% 
better  than  simply  using  the  average  value  of  the  observation.  In  validating 
watershed  prediction  models  that  have  been  calibrated,  Arnold  et  al.  (1997)  term 
values  of  E  greater  than  0.70  as  “a  reasonable  fit,”  while  Micovic  and  Quick 
relate  that  values  of  E  greater  than  0.80  are  “quite  satisfactory.”  For  models  that 
do  not  require  major  calibration,  Frankenberger  et  al.  (1999)  considered  values  of 
E  near  0.60  to  be  “good”. 

While  it  is  recognized  that  no  one  calculation  is  capable  of  quantifying  model 
goodness-of-fit,  a  combination  of  the  above  statistics  should  be  useful.  Legates 
and  McCabe  (1999)  suggest  that  a  statistical  assessment  of  model  performance 
should  include  both  absolute  measures  of  error  (e.g.  B,  RMSE)  and  relative 
measures  of  error  (e.g.  E,  B'2).  Previous  use  of  statistics  in  the  context  of 
reactive  transport  model  validation  will  be  discussed  in  the  following  section. 
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2.4  Case  Studies  Using  Models  to  Simulate  the  Natural 
Attenuation  of  Chlorinated  Solvents 

Several  recent  field  studies  have  focused  on  monitored  natural  attenuation  as  a 

viable  containment  technology  for  chlorinated  solvent  plumes.  Almost  all  of 

these  studies  have  incorporated  reactive  transport  modeling  in  order  to 

demonstrate  that  natural  attenuation  has  taken  place.  These  modeling  efforts 

typically  involve  some  level  of  model  validation  in  order  to  convince  regulators, 

scientists,  and  the  general  public  that  the  application  of  the  model  was  sound. 

This  section  will  review  field  studies  in  the  literature  that  include  model 

calibration,  quantitative  comparative  analysis,  or  qualitative  comparative  analysis 

as  a  means  of  model  validation. 

Model  calibration  is  the  most  oft-performed  validation  technique.  One  of  the  best 
examples  of  a  reactive  transport  model  calibration  can  be  found  in  the  study 
performed  at  the  Area-6  site  on  Dover  Air  Force  Base,  Delaware  (Clement  et  al., 
2000).  RT3D  was  used  in  conjunction  with  the  groundwater  flow  code 
MODFLOW  to  simulate  chlorinated  ethene  fate  and  transport.  Hydrogeologic 
parameter  values  were  determined  through  tracer  tests  or  from  relevant 
literature,  and  contaminant  source  loading  was  quantified  through  model 
calibration.  Model  calibration  was  performed  in  a  trial-and-error  process  to  fit  the 
concentration  profiles  observed  in  1997  for  PCE,  TCE,  cis-DCE,  VC,  ethene,  and 
chloride.  Similar  calibration  studies  include  those  performed  by  Swanson  (1999) 
at  Site  CCFTA-2  at  Cape  Canaveral,  Florida  using  MODFLOW  and  MT3D; 
Moutoux  and  Hicks  (1999)  at  Building  301  at  Offutt  Air  Force  Base,  Nebraska 
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using  MODFLOW  and  MT3D;  Carey  et  al.  (1999)  at  Site  FT-002  at  Plattsburgh 
Air  Force  Base,  New  York  using  MODFLOW  and  BIOREDOX;  Barton  et  al. 
(2000)  at  site  FTP  at  Naval  Air  Station  Fallon,  Nevada  using  MODFLOW  and 
RT3D;  Mason  et  al.  (2000)  at  Naval  Air  Engineering  Station  Lakehurst,  New 
Jersey  using  MODFLOW  and  RT3D. 

In  many  of  the  previously  discussed  studies,  the  calibrated  model  was  used  to 
predict  future  contaminant  concentration  profiles  in  order  to  support  decisions  on 
the  use  of  monitored  natural  attenuation.  The  value  of  these  predictions, 
however,  is  uncertain,  as  predictions  were  not  tested  against  observations. 
Although  the  models  were  calibrated,  there  was  no  assurance  that  the  models 
would  perform  adequately  outside  the  calibration  data  set.  Indeed,  a  study  of 
groundwater  flow  model  calibrations  revealed  that  “good  calibration  does  not  lead 
to  good  prediction”  (Freyberg,  1988).  This  reemphasizes  the  importance  of 
performing  comparative  analysis  between  predictions  to  observations.  While 
examples  can  be  found  in  the  realm  of  weather  forecasting  (Wilks,  1995), 
pesticide  leaching  modeling  (Bond,  1998),  and  groundwater  flow  modeling  (Lee 
and  Ketelle,  1988;  Anderson,  1992;  Bobba,  1993;  Osei,  1995;  Eggleston  and 
Rojstaczer,  2000),  no  documented  studies  involving  chlorinated  ethene  fate  and 
transport  modeling  were  found  that  performed  comparative  analysis  between 
simulations  and  observations  beyond  model  calibration.  If  such  an  analysis  were 
to  be  accomplished,  sufficient  data  would  be  required  to  support  both  model 
calibration  and  comparative  analysis.  Field  sites  that  could  support  such  a 
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demand  for  data  include  Site  FTP  at  Naval  Air  Station  Fallon,  Nevada  (Barton  et 
al.,  2000)  and  Site  LF-01  at  Warren  Air  Force  Base,  Wyoming  (Parsons 
Engineering  Science,  1999). 
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3.0  Methodology 


3.1  Overview 

In  this  chapter,  we  present  a  methodology  for  validating  a  model  that  simulates 
natural  attenuation  of  chlorinated  ethenes.  To  demonstrate  this  methodology,  we 
will  first  choose  a  fate  and  transport  code  that  simulates  the  important  aquifer 
processes  thought  to  be  occurring  at  Site  LF-03  at  F.  E.  Warren  Air  Force  Base, 
Cheyenne,  Wyoming.  Second,  the  code  will  be  calibrated  to  site  data,  and  the 
calibrated  fit  will  be  evaluated  through  the  inspection  of  concentration  contour 
plots  and  goodness-of-fit  statistics.  Third,  the  calibrated  site  model  will  be  used 
to  perform  a  predictive  simulation.  Lastly,  the  simulation  data  will  be  compared 
to  the  data  from  Site  LF-03  using  the  same  evaluation  tools  (e.g.  goodness-of-fit 
statistics)  employed  to  evaluate  calibration.  It  is  hoped  that  this  methodology  will 
increase  confidence  in  our  ability  to  simulate  natural  attenuation  of  chlorinated 
ethenes,  as  well  as  provide  a  greater  understanding  of  natural  attenuation 
processes  occurring  at  a  particular  site. 

3.2  Code  Selection 

3.2.1  Characteristics  of  the  Site 

In  order  to  select  an  appropriate  computer  code  to  simulate  conditions  at  Site  LF- 
03,  it  is  necessary  to  have  a  conceptual  model  of  the  site.  In  this  section,  we 
briefly  characterize  the  site,  to  include  a  description  of  the  military  base,  the 
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landfill  under  investigation,  and  the  subsurface  processes  thought  to  be  taking 
place  near  that  landfill. 

F.  E.  Warren  Air  Force  Base  (Warren  AFB),  adjacent  to  the  city  limits  of 
Cheyenne,  Wyoming,  has  been  a  military  post  since  1867.  The  Air  Force 
assumed  control  of  the  installation  in  1 947,  and  has  used  the  base  to  support 
both  training  and  operational  commands.  Warren  AFB  is  currently  the  home  of 
Air  Force  Space  Command’s  90th  Space  Wing,  a  unit  responsible  for  the 
readiness  and  maintenance  of  intercontinental  ballistic  missiles  (ICBMs). 

Landfill  03  (LF-03)  was  operational  on  Warren  AFB  from  approximately  the  mid- 
1950’s  to  the  mid-1960’s  for  the  disposal  of  industrial  and  residential  wastes 
generated  by  the  base  (Parsons  Engineering  Science,  1999).  LF-03  is  located  in 
the  southeast  corner  of  the  base,  is  approximately  7  acres  in  area,  and  contains 
a  maximum  of  15,400,000  cubic  ft  (ft3)  of  fill.  In  the  mid  1980’s  it  was  determined 
that  leachate  from  LF-03  could  possibly  pose  a  threat  to  human  health  and  the 
environment.  Groundwater  monitoring  wells  were  installed  and  sampled  in  1987 
and  1988.  After  Warren  AFB  was  placed  on  the  National  Priorities  List  (NPL)  in 
February  of  1990,  a  1991  remedial  investigation  (Rl)  revealed  that  TCE  was  the 
primary  groundwater  contaminant.  It  was  hypothesized  that  TCE,  used  as  a 
metal  cleaner  at  base  shops,  had  been  disposed  in  LF-03  according  to  the 
standard  industry  practices  of  the  day.  A  focused  Rl  was  started  for  LF-03  in 
1995,  and  a  treatability  study  to  evaluate  the  use  of  monitored  natural  attenuation 
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(MNA)  was  conducted  by  Parsons  Engineering  Science  and  the  USEPA  National 
Risk  Management  Research  Laboratory  (NRMRL)  in  1999. 

The  treatability  study  indicated  that  several  processes  were  taking  place  in  the 
aquifer.  First,  it  suggested  that  anaerobic  conditions  near  the  source  area 
(caused  by  the  degradation  of  co-contaminants)  allowed  for  the  reductive 
dechlorination  of  TCE.  Several  lines  of  evidence  supported  this  claim,  including 
the  presence  of  cis-1 ,2-DCE.  As  DCE  was  not  used  in  the  base  shops,  the  most 
likely  source  of  this  contaminant  is  from  the  reduction  of  TCE.  Second,  it  was 
suggested  that  DCE  and  VC  were  being  oxidized  downgradient  of  the  source. 

The  most  convincing  evidence  for  this  theory  was  the  increasing  ratio  of  TCE  to 
DCE  found  downgradient  of  the  landfill.  This  led  to  the  hypothesis  that  “DCE  is 
degraded  through  oxidation  reactions,  while  TCE  mass  is  relatively  unaffected  by 
destructive  attenuation  mechanisms”  (Parsons  Engineering  Science,  1999). 

Oxygen,  present  in  significant  concentrations  downgradient  of  the  source,  is 
assumed  to  be  the  electron  acceptor  in  the  oxidation  reactions.  Unfortunately, 
there  are  insufficient  oxygen  data  to  model  direct  oxidation  of  DCE  and  VC.  As  a 
result,  it  will  be  assumed  that  the  attenuation  of  all  chlorinated  ethenes  is  the 
result  of  first-order  reductive  dehalogenation. 

Assumptions  are  also  made  about  the  sorption  process.  Following  the  treatability 
study  (Parsons  Engineering  Science,  1999),  sorption  is  assumed  to  be  an 
equilibrium  process.  The  linear  sorption  model  is  also  deemed  appropriate  for 
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this  site,  as  the  concentration  of  contaminants  is  relatively  low  throughout  the 
area  under  study.  Based  on  these  assumptions,  any  computer  code  chosen  for 
this  research  should  be  able  to  simulate  linear  equilibrium  sorption  and  reductive 
dehalogenation  of  chlorinated  ethenes. 


3.2.2  Criteria  for  Code  Selection 

When  deciding  which  computer  code  to  apply  to  a  site,  the  most  important 
consideration  is  that  “the  processes  identified  as  being  important  at  the  site  must 
correspond  to  those  included  in  the  model”  (Weaver  et  al.,  1987).  The  computer 
code’s  ability  to  model  a  relevant  process  is  known  as  an  essential  code 
capability,  in  that  “if  a  candidate  code  does  not  include  the  essential  capabilities, 
it  should  be  removed  from  consideration”  (ASTM,  1997).  Due  to  the  complex 
flow  characteristics  found  at  Site  LF-03  (Parsons  Engineering  Science,  1999),  it 
was  decided  that  a  3D  transport  model  was  necessary  to  represent  the  site.  3D 
representation  is  considered  an  essential  code  capability,  as  is  the  ability  to 
represent  linear  equilibrium  sorption  and  anaerobic  reduction.  Table  3.1 
compares  essential  code  capabilities  against  the  capabilities  of  the  candidate 
codes. 


Essential  Code  Capabilities 

Biochlor 

Bio-Redox 

BR3D 

RT3D  | 

IPhvsicochemical  Processes  I 

3-D  Advection 

No 

Yes 

Yes 

Yes 

3-D  Dispersion 

Yes 

Yes 

Yes 

Linear  Equilibrium  Sorption 

Yes 

Yes 

Yes 

Bioloaical  Processes 

(Anaerobic  Reduction  of  chlorinated  ethenes 

Yes  Yes 

Yes 

Yes  | 

Table  3.1  Essential  Capabilities  of  Candidate  Codes 
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Based  upon  the  essential  code  capabilities,  only  Bio-Redox,  BR3D,  and  RT3D 
will  be  considered  for  further  study. 


Model  selection  also  depends  upon  non-essential  code  capabilities,  such  as 
ease  of  use  and  code  credibility.  These  capabilities  are  subjectively  ranked,  in 
order  of  decreasing  relative  importance,  by  the  modeler.  This  ranking  process 
allows  a  balance  to  be  struck  between  modeling  effort  and  results  obtained 
(ASTM,  1997).  Table  3.2  contains  a  list  of  preferred  attributes  on  which  the 
candidate  codes  can  be  judged.  The  ability  of  each  program  to  meet  designated 
non-essential  capability  is  ranked  Good,  Fair,  or  Poor.  User  Support  and  Ease  of 
Use  are  relatively  important  in  this  effort  due  to  the  experience  of  the  modeler 
and  the  time  available  in  the  study.  As  this  research  will  attempt  to  build 
confidence  in  the  code’s  application  to  the  site,  attributes  such  as  Code 
Acceptance  and  Code  Credibility  are  less  important. 


Non-Essential  Code  Capabilities 

Bio-Redox 

BR3D 

RT3D 

User  Support 

Poor 

Good 

Poor 

Ease  of  Use  (data  input/output) 

Poor 

Poor 

Poor 

Code  Documentation 

Poor 

Poor 

Good 

Code  Availability 

Good 

Good 

Code  Acceptance 

Good 

Poor 

Good 

Code  Credibility 

Good 

Poor 

Good 

Table  3.2  Non-Essential  Capabilities  of  Candidate  Codes 


3.2.3  Code  Selected  to  Represent  Site  LF-03 

Of  the  codes  inventoried  in  Chapter  2,  BIOREDOX,  RT3D,  and  BR3D  are  able  to 
perform  the  essential  calculations  required  to  represent  the  physicochemical  and 
biological  processes  thought  to  be  occurring  at  Site  LF-03.  This  results  in 
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choosing  a  computer  code  based  on  the  non-essential  capabilities  listed  in  the 
previous  paragraph.  Of  these  capabilities,  the  most  important  is  the  availability  of 
technical  support  to  the  modeler.  Based  on  these  criteria,  it  was  decided  to  use 
BR3D  in  the  validation  studies  performed  in  this  research. 

3.3  Code  Application  to  Site  LF-03 

The  application  of  BR3D  to  Site  LF-03  is  based  upon  the  treatability  study  (TS)  of 
MNA  performed  by  Parsons  Engineering  Science  (1999).  In  the  TS,  the 
computer  code  MT3D  was  used  in  conjunction  with  MODFLOW  (a  groundwater 
flow  model)  to  depict  the  fate  and  transport  of  TCE.  A  220-cell  by  125-cell  by 
three  layer  model  domain  was  created  (see  Figure  3.1),  with  each  cell  measuring 
20  feet  by  20  feet.  MODFLOW  and  MT3D  were  calibrated  using  data  collected 
from  27  monitoring  wells.  Simulations  were  then  performed  to  ascertain  the 
effectiveness  of  different  remedial  alternatives.  While  this  modeling  study  proved 
informative  to  decision-makers  at  Warren  AFB,  it  had  several  shortcomings. 

First,  TCE  was  the  only  contaminant  studied,  even  though  DCE  and  VC  were 
also  detected  at  the  site.  Second,  the  study  did  not  involve  any  effort,  outside  of 
model  calibration,  to  validate  the  application  of  the  model  to  the  site.  This 
research  will  address  these  two  points  specifically  by  using  BR3D  to  predict  TCE, 
DCE,  and  VC  concentrations,  then  compare  those  values  to  values  observed  at 
the  site.  The  remainder  of  this  section  will  discuss  how  the  groundwater  flow 
model  created  in  the  TS  was  adapted  for  use  in  this  research,  and  will  describe 
how  the  contaminant  fate  and  transport  code  BR3D  was  applied  to  LF-03. 
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Figure  3.1  Model  Domain  (after  Parsons  Engineering  Science  (1999)) 


3.3.1  Groundwater  Flow  Model 

Like  MT3D,  BR3D  requires  the  use  of  MODFLOW,  a  groundwater  flow  code 
developed  by  the  U.S.  Geological  Survey  (Harbaugh  and  McDonald,  1996). 
MODFLOW  calculates  hydraulic  heads  and  groundwater  fluxes,  which  are  then 
used  by  fate  and  transport  codes  such  as  MT3D  and  BR3D  to  calculate 
contaminant  concentrations.  MODFLOW  is  a  highly  regarded  program,  and  its 
extensive  use  in  the  study  of  hydrogeology  and  contaminant  transport  bears 
witness  to  its  usefulness.  This  research  does  not  attempt  to  validate  the 
application  of  MODFLOW  to  Site  LF-03.  Instead,  the  Visual  MODFLOW  (Version 
2. 8.2.0)  flow  model  of  Site  LF-03  (Parsons  Engineering  Science,  1999)  is 
assumed  to  be  valid  for  the  purposes  of  this  research.  The  site  flow  model  is 
based  upon  reasonable  assumptions  and  a  thorough  conceptual  understanding 
of  the  site  (Parsons  Engineering  Science,  1999).  Figure  3.2  depicts  hydraulic 
head  contours  at  the  site  based  on  measurements  taken  May  1999.  The 
reported  hydraulic  mass  balance  bolsters  the  assumption  of  flow  model  validity; 
the  discrepancy  between  incoming  and  outgoing  hydraulic  flux  for  the  steady- 
state  calibrated  flow  model  was  calculated  as  0.0  percent  (Parsons  Engineering 
Science,  1 999).  Groundwater  flow  parameter  values  that  apply  to  the  entire 
model  domain  are  listed  in  Table  3.3,  while  individual  cell  values  can  be  found  in 
the  flow  model  itself  (Parsons  Engineering  Science,  1999).  The  average  value  of 
hydraulic  conductivity  depends  upon  which  of  the  three  zones  (ungradient, 
downgradient,  or  Crow  Creek)  is  under  consideration  (Table  3.3  and  Figure  3.2). 
As  it  is  assumed  that  the  groundwater  flow  portion  of  the  model  will  perform 
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CROW  CREEK  ZONE 


Figure  3.2  Shallow  Groundwater  Surface  (after  Parsons  Engineering  Science  (1999)) 


Variable 

Value 

Estimated  Effective  Porosity 

0.2 

Bulk  Density  of  Aquifer  Solids 

1 .65  kg/I 

Average  Hydraulic  Gradient 

0.02  ft/ft 

Average  Hydraulic  Conductivity 

Range 

Crow  Creek  Zone  (see  Figure  3.2) 

Average  Groundwater  Velocity 

iim— WM— 

Range 

Geometric  Average 

6.9  ft/yr  | 

Table  3.3  Groundwater  Flow  Parameter  Values 

adequately,  this  research  will  focus  on  the  application  of  the  contaminant 
transport  model  to  the  site. 


3.3.2  Contaminant  Fate  and  Transport  Model 

This  research  will  model  the  natural  attenuation  TCE,  DCE,  and  VC  due  to 
reductive  dehalogenation.  This  section  will  discuss  how  BR3D  will  be  applied  to 
the  site,  to  include  fate  and  transport  processes  modeled  and  initial  and 
boundary  conditions  assumed. 


3.3.2. 1  Processes  Modeled 

BR3D  can  model  3D  advection,  3D  dispersion,  equilibrium  or  rate-limited 
sorption,  and  several  biodegradation  processes.  Assuming  steady  state 
groundwater  flow,  linear  sorption,  and  constant  porosity  (both  spatially  and 
temporally),  these  processes  can  be  incorporated  into  the  general  transport 
equation  for  a  single  contaminant  (Charbeneau,  2000;  Domenico  and  Schwartz, 
1998): 
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(3.1) 


R^2  =  _v.vC  +  V.(D.VC)  +  rbi0 
dt 

where 

C  =  concentration  of  contaminant  in  aqueous  phase  [M/L3] 
t  =  time  [T] 

v  =  average  linear  velocity  vector  (in  the  x,y,  and  z  direction)  [L/T] 

D  =  dispersion  coefficient  matrix  [L2/T] 
rbi0  =  source/sink  term  for  contaminant  production/destruction  in  aqueous  phase 
[M/L3-T] 

R  =  retardation  factor  of  contaminant  due  to  sorption  [-] 

On  the  right  hand  side  of  the  equation,  the  first  term  represents  advection,  the 
second  dispersion,  and  the  third  biological  reactions.  Linear  equilibrium  sorption 
is  included  on  the  left  side  of  the  equation  in  the  R  term.  The  remainder  of  this 
section  will  describe  the  specific  sorption  and  reaction  equations  used  in  BR3D  to 
solve  the  general  transport  equation. 

Retardation  due  to  linear  equilibrium  sorption  can  be  described  using  the 
equation 

R=1+Kd^  (3.2) 

where 

R  =  retardation  factor  of  contaminant  [-] 

Kd  =  distribution  coefficient  of  contaminant  (L3/M) 
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0  =  porosity  (-) 

pb  =  bulk  density  of  aquifer  solids  (M/L3) 

The  distribution  coefficient  of  the  contaminant,  KD  ,  can  be  determined  by 
(Domenico  and  Schwartz,  1998): 

=  Kocfoc  (3-3) 

where 

Koc  =  partition  coefficient  of  a  contaminant  between  organic  carbon  and 
water  (L3/M) 

foc  =  weight  fraction  of  organic  carbon  in  the  aquifer  solids  (-) 

As  Koc  values  were  not  experimentally  determined  at  Site  LF-03,  they  may  be 
estimated  using  the  relationship  (Domenico  and  Schwartz,  1998) 

l°gK0C  =-0.21  +  logKow  (3.4) 

where 

Kow  =  partition  coefficient  of  a  contaminant  between  octanol  and  water  (-) 

BR3D  also  models  the  biodegradation  processes  at  Site  LF-03.  Reductive 
dehalogenation  of  chlorinated  ethenes  is  assumed  to  be  the  dominant  process, 
and  is  described  using  first  order  kinetics: 

rli0=^=-kC  (3.5) 

where 

C  concentration  of  dissolved  contaminant  (M/L3) 
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k  =  contaminant  decay  first  order  rate  constant  (T'1) 


3.3. 2.2  Initial  and  Boundary  Conditions 

It  is  assumed  that  TCE  started  leaching  into  the  groundwater  in  1960,  and  that  no 
contamination  was  present  before  that  time.  TCE  is  the  only  contaminant 
infiltrating  into  the  groundwater;  DCE  and  VC  are  assumed  to  be  present  from 
the  reduction  of  aqueous  phase  TCE.  It  is  also  assumed  that  TCE  infiltrates  from 
cells  in  three  designated  source  zones,  as  shown  in  Figure  3.3.  The  mass 
loading  rate  of  TCE  entering  the  aquifer  is  listed  in  Table  3.4. 


Time 

Source  Zone  1 

Source  Zone  2 

Source  Zone  3 

1960-1965 

0.0028 

0 

0 

1965-1970 

0.0028 

0 

0 

1970-1975 

0.0028 

0 

0 

1975-1980 

0.0019 

0.0094 

0 

1980-1985 

HKESIEHH 

0.0047 

0.0094 

1985-1990 

ITiU 

0.0024 

0 

1990-1995 

0.0006 

0.0094 

0 

1995-1999 

0.0004 

0.0047 

0 

Table  3.4  TCE  Mass  Loading  per  Cell  (in  kg/yr)  (after  Parsons  Engineering 

Science,  1999) 


3.4  Site  Model  Validation 

In  this  research,  model  calibration  will  be  performed  using  data  collected  in  1993, 
while  a  predictive  simulation  will  be  compared  to  observations  made  in  1999. 
Both  calibration  and  validation  will  include  concentration  contour  plots  and 
goodness-of-fit  statistics  in  order  to  evaluate  model  performance. 
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208  _rk_  MONITORING  WELL  LOCATION 


Figure  3.3  Model  Source  Area  Configuration  (after  Parsons  Engineering  Science  (1999)) 


3.4.1  Calibration 


BR3D  will  be  calibrated  using  concentrations  of  TCE,  DCE,  and  VC  observed  in 
1993  (see  Table  3.5)  through  the  use  of  RMSE. 


Site 

Well# 

Corresponding  Model  Location 

Obs  Pt  Row  Col  Layer  X  Y 

Observed  Concentrations 

TCE  cis-1 ,2DCE  VC 

62 

1 

46 

87 

1 

1736 

1588 

2.00E-05 

1.22E-05 

0.00E+00 

147 

2 

104 

1 

2065 

1611 

1.50E-05 

8.20E-06 

0.00E+00 

198 

3 

32 

103 

1 

2045 

1866 

0.00E+00 

199 

4 

48 

1 

EHH 

1553 

1 .20E-06 

O.00E+00 

201 

6 

69 

1 

1139 

1.30E-06 

O.OOE+OO 

203 

7 

70 

119 

1 

2375 

1108 

8.80E-06 

O.OOE+OO 

0.00E+00 1 

206 

8 

90 

1 

1139 

6.40E-06 

1.00E-06 

9 

1 

1365 

1.10E-05 

208 

10 

1 

1436 

1509 

2.00E-05 

1.60E-05 

0.00E+00 

209 

11 

55 

83 

1 

3.40E-05 

2.80E-05 

0.00E+00 

210 

13 

53 

105 

1 

4.40E-05 

2.40E-05 

9.00E-07 

15 

ng 

1 

1 .90E-05 

9.10E-06 

9.00E-07 

16 

■s 

k m 

1 

BcliBfSl 

FE3 

0.00E+00 

O.OOE+OO 

17 

mm 

50 

1 

989 

1393 

0.00E+00 

O.OOE+OO 

0.00E+00 

|  236 

18 

64 

150 

1 

2998 

1227 

0.00E+00 

O.OOE+OO 

0.00E+00 

Table  3.5  Observed  Chlorinated  Ethene  Concentrations  for  1993  (after 
Parsons  Engineering  Science  (1999)) 


By  inspecting  the  values  observed  in  1993,  it  can  be  seen  that  chlorinated  ethene 
concentrations  are  much  larger  at  some  observation  points  than  at  others.  For 
example,  the  concentration  of  TCE  at  Observation  Point  13  is  over  35  times 
greater  than  TCE  concentration  at  Observation  Point  4.  If  we  fit  the  raw  data, 
these  large  values  will  exert  an  undue  influence  on  the  calibration  by  masking  the 
contribution  of  other  observation  points  to  the  overall  model  fit.  In  order  to 
dampen  the  effect  of  these  extreme  values  on  the  goodness-of-fit  statistics,  it 
was  decided  to  perform  a  log  transformation  on  the  concentration  values.  As  the 
data  values  themselves  are  of  interest  (as  opposed  to  the  difference  in  values 
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being  important)  and  due  to  the  large  range  of  the  values,  logs  of  base  10  will  be 
used  (Cole,  2000).  The  transformed  data  will  then  be  analyzed  using  RMSE. 


The  parameters  to  be  varied  in  the  calibration  are  listed  in  Table  3.6,  while 
parameters  held  constant  are  listed  in  Table  3.7. 


Parameter 

Symbol 

Ranae  of  Values 

Source 

First  order  rate  constant  for  TCE 

Ktce 

0  -  0.023  d'1 

Suarez  and  Rifai,  1999 

First  order  rate  constant  for  DCE 

1<DCE 

0- 0.130  d'1 

Suarez  and  Rifai,  1999 

First  order  rate  constant  for  VC 

kvc 

0  -  0.007  d’1 

Suarez  and  Rifai,  1999 

Table  3.6  Model  Parameters  Varied  During  Calibration 


Parameter 

Svmbol 

Value 

Source 

Longitudinal  Dispersivity 

aL 

25  ft 

Parsons  Engineering  Science,  1999 

Transverse  Dispersivity 

aT 

2.5  ft 

Parsons  Engineering  Science,  1999 

Vertical  Dispersivity 

av 

0.9  ft 

Parsons  Engineering  Science,  1999 

Retardation  factor  for  TCE 

Rtce 

1.33 

CRC  Press,  1997;  Equations  3.2  -  3.4 

Retardation  factor  for  DCE 

1.28 

CRC  Press,  1997;  Equations  3.2  -  3.4 

Retardation  factor  for  VC 

CRC  Press,  1997;  Equations  3.2  -  3.4 

Table  3.7  Model  Parameters  Held  Constant  During  Calibration 

The  first  order  decay  rate  constants  for  TCE,  DCE,  and  VC  will  be  varied  to 
minimize  the  RMSE  of  all  three  contaminants.  Observed  and  calibrated 
contaminant  contour  plots  and  optimized  RMSE,  bias(B)  and  coefficient  of 
efficiency  (E)  will  be  reported  for  each  contaminant  as  a  measure  of  model 
calibration. 
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3.4.2  Quantitative  Comparative  Analysis 

After  BR3D  is  calibrated  to  1993  data,  a  predictive  simulation  will  be  performed. 
The  model  will  be  run  from  1960  to  1999,  and  results  will  be  compared  to  the 
1999  observed  concentrations  of  TCE,  DCE,  and  VC.  As  with  the  model 
calibration,  observed  and  calibrated  contaminant  contour  plots  and  RMSE,  bias 
(B)  and  coefficient  of  efficiency  (E)  will  be  reported  for  each  contaminant. 
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4.0  Analysis 


4.1  Calibration  Results 

Chlorinated  ethene  decay  rate  constants  were  estimated  using  a  forward,  trial- 
and-error  process.  Using  the  sum  of  the  three  chlorinated  ethene  RMSEs  as  the 
objective  function,  the  decay  rates  for  DCE  and  VC  were  held  constant  while  the 
decay  rate  for  TCE  was  varied  to  minimize  the  RMSE  sum  (see  Figure  4.1). 


Figure  4.1  Calibration  of  TCE  First  Order  Decay  Constant 
for  kDCE=1x10-5  d'1  ,  kvc  =  1X10-6  d'1 

Once  a  kTCEwas  found  that  minimized  the  objective  function  (sum  of  RMSE), 
kDcE  was  varied  while  kTcE  and  kVc  were  held  constant.  This  process  was 
repeated  until  kTcE,  kDCE,  and  kVc  each  converged  to  a  value  that  produced  the 
(hopefully)  global  minimum  objective  function.  The  calibrated  decay  rates  are 
listed  in  Table  4.1. 
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First  Order  Rate  Constant,  k 

(4-1) 

2.3x10-4 

1.8x10-4 

VC 

0 

Table  4.1  Calibrated  Model  Parameters 

The  fit  values  for  first  order  decay  are  reasonable,  in  that  they  fall  within  the 
range  of  values  reported  in  the  literature  (Suarez  and  Rifai,  1999).  It  can  also  be 
seen  that  the  greater  number  of  chlorine  substituents  a  contaminant  has,  the 
greater  the  rate  of  decay  (kTcE>kDCE>kvc)-  These  two  facts  support  the 
assumption  that  the  aquifer  is  anaerobic,  and  that  reductive  dehalogenation  is  an 
important  process  at  the  site. 

Significant  insight  can  be  gained  through  qualitative  evaluation  of  the  calibrated 
fit.  To  this  end,  simulated  and  observed  concentration  contour  plots  for  TCE 
(Figures  4.2  and  4.3,  respectively),  DCE  (Figures  4.4  and  4.5),  and  VC  (Figures 
4.6  and  4.7)  were  generated  using  the  linear  interpolation  option  of  Surfer 
(Version  5.03),  a  surface  mapping  program. 

In  comparing  the  concentration  contour  plots,  we  found  that  simulated 
concentrations  matched  observed  concentrations  reasonably  well  for  TCE  and 
VC,  while  the  DCE  concentration  data  were  less  well  matched.  Specifically,  it  is 
apparent  that  the  calibrated  model  slightly  underestimates  the  observed 
concentration  of  DCE  found  at  Site  LF-03. 
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Figure  4.2  Calibrated  TCE  Concentration  Contour  Plot,  Layer  1  (g/L) 
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Figure  4.4  Calibrated  DCE  Concentration  Contour  Plot,  Layer  1  (g/L) 


Figure  4.5  1993  Observed  DCE  Concentration  Contour  Plot,  Layer  1  (g/L) 
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Figure  4.6  Calibrated  VC  Concentration  Contour  Plot,  Layer  1  (g/L) 
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Figure  4.7  1993  Observed  VC  Concentration  Contour  Plot,  Layer  1  (g/L) 
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Goodness-of-fit  statistics  were  also  used  to  evaluate  the  calibrated  fit.  Table  4.2 
contains  statistics  that  describe  the  absolute  magnitude  of  errors  (RMSE), 
relative  magnitude  of  errors  (E),  and  direction  and  magnitude  of  the  model’s 
tendency  for  over-  or  under-prediction  (B,  B'2)  for  each  chlorinated  ethene. 


RMSE 

mm 

mm 

(log(ppm)) 

■H 

BH 

TCE 

0.71 

-0.27 

DCE 

0.74 

0.51 

-0.42 

0.16 

VC 

0.21 

0.60 

-0.030 

0.008 

Table  4.2  Goodness-of-Fit  Statistics  for  Model  Calibration 


The  model  bias  (B)  of  each  contaminant  reveals  that  concentrations  for  all  three 
are  being  under-estimated,  with  DCE  having  the  greatest  magnitude  of  bias  (B'2). 
This  under-estimation  could  have  significant  implications  if  the  model  is  to  be 
used  to  predict  chlorinated  ethene  fate  and  transport.  However,  the  coefficient  of 
efficiency  for  the  three  contaminants  (E),  while  less  than  the  ideal  value  of  1 ,  is 
greater  than  zero.  Using  E  as  a  gauge  of  model  performance  is  relatively  new  to 
contaminant  transport  modeling.  Accordingly,  there  is  no  consensus  of  what 
value  of  E  indicates  an  acceptable  fate  and  transport  model,  especially  when  the 
data  have  been  log  transformed.  However,  as  this  model  did  not  require  a  major 
calibration  (only  three  parameters  varied),  it  may  be  possible  to  follow  the 
example  given  by  Frankenberger  et  al.  (1999)  for  the  validation  of  a  hydrology 
model.  If  this  example  is  applicable,  values  in  the  range  of  E  =  0.60  can  be 
considered  good,  as  is  the  case  for  TCE  and  DCE  fit.  Values  of  E  <  0.5,  as  in  the 
case  for  VC  fit,  tend  to  indicate  poor  model  performance.  However,  it  should  be 
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noted  that  the  VC  fit  was  based  on  only  two  quantifiable  data  points  (wells  210 
and  21 1  in  Table  3.5),  and  that  the  kVc  is  essentially  zero.  This  essentially 
means  that  only  two  fitting  parameters  (kTcE  and  kDcE)  were  used  to  calibrate  the 
model  to  TCE  and  DCE  data.  This  fact  was  taken  into  consideration  along  with 
the  evaluation  of  concentration  contour  plots  and  the  goodness-of-fit  statistics. 
After  a  review  of  the  information,  we  concluded  that  the  calibrated  model 
provided  a  good  fit  to  the  TCE  and  DCE  data,  and  an  acceptable  fit  to  the  VC 
data.  The  fact  that  the  model  can  be  reasonably  calibrated  with  minimal  fitting 
parameters  suggests  that  the  modeling  assumptions  of  linear  equilibrium  sorption 
and  first  order  decay  are  not  invalid. 

4.2  Predictive  Simulation  and  Comparative  Analysis 

After  the  decay  rate  constants  were  calibrated,  a  predictive  simulation  was  run 
from  1960  to  1999.  Simulation  values  were  then  compared  to  values  observed  at 
Site  LF-03  in  1999  (Table  4.3).  Figures  4.8,  4.9,  and  4.10  show  the  simulated 
and  observed  chlorinated  ethene  concentrations. 

Inspection  of  the  4x1  O'5  contour  lines  in  Figure  4.8  indicates  that  simulated  TCE 
concentration  matches  well  to  observed  data,  while  an  inspection  of  Figure  4.9 
shows  that  DCE  is  considerably  underestimated  in  the  simulation.  In  contrast, 
the  model  overestimates  VC  concentrations  (see  Figure  4.10).  It  should  be 
noted,  however,  that  the  simulated  VC  concentration  are  nearly  equal  to  the  VC 
detectable  limit  of  1x1  O'7  g/L,  such  that  the  overestimation  is  relatively  small. 
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Site 

Well# 


62 


47 

98 
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199M1 

201 


203 

206 


_2 

2 


m 


Corresponding 

Obs  Pt  Row  Col 


46  87 


45 


48  161 


48 

69 


119 


90 


57  75 

50  72 


55  83 

55  83 


Model  Location 

Layer  X  Y 


1  1736  1588 


1611 


1866 


3210  1553 


553 


inrr 


64  150 

81  113 


77  133 


53  74 


PES-2D 

PES-3D 


0  1450 
0  1450 


2349  1500 
3088  1786 
989  1393 


1227 

895 


2655 

1473 


759  |  1591 
786 


Observed  Concentrations  (g/L) 

TCE 

Cis-1,2DCE 

VC 

5.00E-07 

5.00E-07 

0.00E+00 

1.04E-05 

5.90E-06 

O.00E+00 

0.00E+00 

0.00E+00 

0.00E+00 

0.00E+00 

0.00  E+00 

0.00E+00 

0.00E+00 

0.00E+00 

0.00E+00 

5.00E-07 

0.00E+00 

0.00E+00 

1.28E-05 

5.10E-06 

0.00E+00 

2.60E-06 

1.20E-06 

0.00E+00 

EBEEa 

3.66E-05 

0.00E+00 

2.28E-05 

0.00E+00 

9.31  E-05 

1.29E-04 

7.71  E-05 

5.1  IE-05 

O.00E+00 

UsVJJilil 

1.59E-05 

0.00E+00 

IRWaSVIdl 

6.90E-06 

O.OOE+OO 

3.30E-06 

0.00E+00 

0.00E+00 

0.00E+00 

0.00E+00 

5.00E-07 

0.00E+00 

0.00E+00 

5.00E-07 

0.00E+00 

0.00E+00 

0.00  E+00 

0.00E+00 

0.00E+00 

1.12E-06 

0.00E+00 

0.00E+00 

3.39E-05 

3.35E-05 

0.00E+00 

4.30E-06 

2.30E-06 

0.00E+00 

5.00E-07 

0.00E+00 

0.00E+00 

1.23  E-05 

3.00E-06 

0.00E+00 

3.30E-06 

1.30E-06 

0.00E+00 

1.10E-05 

0.00E+00 

Table  4.3  Observed  Chlorinated  Ethene  Concentrations  for  1999  (after 
Parsons  Engineering  Science  (1999)) 
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Figure  4.8  Simulated  (a)  and  Observed(b)  TCE  Concentrations  for  1999 


Figure  4.9  Simulated  (a)  and  Observed  (b)  DCE  Concentrations  for  1999 


X  Coordinate  (ft)  X  Coordinate  (ft) 

Figure  4.10  Simulated  (a)  and  Observed  (b)  VC  Concentrations  for  1999 
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To  complete  the  evaluation  of  model  performance,  goodness-of-fit  statistics  for 
the  predictive  simulation  (summarized  in  Table  4.4)  were  calculated. 


RMSE 

(log(ppm)) 

mm 

B 

(log(ppm)) 

mm 

1  TCE 

0.50 

0.77 

-0.09 

0.009 

0.77 

0.47 

-0.39 

0.14 

0.24  | 

- 

0.13 

- 

Table  4.4  Goodness-of-Fit  Statistics  for  1999  Data 


As  VC  was  not  observed  at  Site  LF-03  in  1999,  the  observed  variance  of  VC 
concentrations  equals  zero.  This  prevents  the  calculation  of  the  Nash-Sutcliffe 
coefficient  of  efficiency  (E)  and  the  non-dimensional  bias  (B'2)  for  this  data  set. 
Regardless,  interesting  observations  can  be  made  from  the  statistics.  First,  the 
magnitude  of  non-dimensional  bias  (B/2)  for  TCE  decreased  from  0.069  in  1993 
to  0.009  in  1999.  It  is  hypothesized  that  the  TCE  loading  rate  in  the  first  few  time 
steps  (see  Table  3.4)  may  have  been  insufficient,  resulting  in  the  underestimation 
of  TCE  concentration  later  in  the  simulation.  This  reduction  in  the  magnitude  of 
bias  has  implications  for  another  goodness-of-fit  statistic,  as  E  is  related  to  B'2  by 
equation  2.18: 

E  =  R2-C2-B'2  (2.18) 

The  decrease  in  B/2  could  account  for  the  unusual  circumstance  that  the  model 
fits  the  observed  data  better  when  in  the  predictive  mode  than  when  the  fit  was 
calibrated,  as  shown  by  the  increase  in  E  from  0.71  to  0.77.  The  value  of  E  for 
TCE  suggests  that  the  TCE  concentration  fit  is  good  (Frankenberger  et  al., 

1999).  Second,  the  value  of  E  for  DCE  is  relatively  low  (below  0.50),  indicating 
that  the  fit  of  simulated  DCE  to  the  observed  values  is  relatively  poor 
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(Frankenberger  et  al.,  1999).  Lastly,  the  positive  value  of  VC  bias  (B)  indicates 
that  VC  is  being  over-predicted.  As  no  VC  was  observed  at  Site  LF-03  in  1999, 
any  model  prediction  over  1x1  O'7  g/L  (the  detection  limit  for  VC)  would  result  in 
over-prediction. 

While  chlorinated  ethene  advection  is  relatively  slow,  and  only  six  additional 
years  have  been  simulated,  we  found  that  BR3D  provided  good  predictions  of 
TCE  concentrations  and  fair  predictions  of  DCE  concentration,  while  lack  of 
quantifiable  data  precludes  the  evaluation  of  VC  prediction.  From  these 
observations  of  model  performance,  we  inferred  that  the  assumptions  of  linear 
equilibrium  sorption  and  reductive  dehalogenation  were  not  invalid  when 
describing  TCE  and  DCE  fate  and  transport. 
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5.0  Conclusion 


5.1  Summary 

In  this  thesis,  the  contaminant  fate  and  transport  model  BR3D  was  applied  to  Site 
LF-03  on  F.E.  Warren  Air  Force  Base  to  simulate  the  natural  attenuation  of  TCE, 
DCE,  and  VC  through  reductive  dehalogenation.  Reductive  dehalogenation  was 
modeled  as  a  first  order  process,  with  decay  rate  constants  being  calibrated 
through  a  forward,  trial-and-error  process.  After  model  calibration,  a  predictive 
simulation  was  performed.  Simulated  values  for  chlorinated  ethene 
concentrations  were  compared  to  observed  values  using  concentration  contour 
plots  and  goodness-of-fit  statistics. 

5.2  Calibration 

The  calibrated  values  for  first  order  decay  were  found  to  be  reasonable,  in  that 
they  fell  within  the  range  of  values  reported  in  the  literature  (Suarez  and  Rifai, 
1999).  Contaminants  with  greater  number  of  chlorine  substituents  had  greater 
rates  of  decay  (kTcE>kDCE>kvc).  as  would  be  expected  in  an  anaerobic  aquifer 
where  reductive  dehalogenation  is  taking  place.  After  a  simulation  was  run  using 
the  fit  decay  rates,  a  comparison  of  the  concentration  contour  plots  revealed  that 
simulated  concentrations  matched  observed  concentrations  reasonably  well  for 
TCE  and  VC,  while  the  simulated  DCE  concentration  tended  to  underestimate 
the  observed  data.  It  should  be  noted  that  choosing  other  decay  rates  could 
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have  improved  DCE  fit,  but  this  would  have  resulted  in  the  degeneration  of  TCE 
and  VC  fit  such  that  the  overall  fit  of  the  model  would  have  suffered. 

The  model  bias  for  each  contaminant  revealed  that  TCE,  DCE,  and  VC  were 
being  under-estimated,  with  DCE  having  the  greatest  magnitude  of  bias. 
However,  the  coefficient  of  efficiency  for  the  three  contaminants  (E),  showed  that 
the  model  performed  better  than  the  mean  value  at  estimating  chlorinated  ethene 
concentrations.  Specifically,  values  of  E  were  judged  as  good  for  TCE  and  DCE 
but  poor  for  VC.  As  BR3D  is  minimally  calibrated  (essentially  only  two  fitting 
parameters)  we  concluded  that  the  calibrated  model  provided  a  good  fit  to  the 
TCE  and  DCE  data,  and  an  acceptable  fit  to  the  VC  data.  The  fact  that  the 
model  can  be  reasonably  calibrated  with  minimal  fitting  parameters  suggests  that 
the  modeling  assumptions  of  linear  equilibrium  sorption  and  first  order  decay  are 
not  invalid. 

5.3  Predictive  Simulation  and  Comparative  Analysis 

After  BR3D  was  calibrated,  a  predictive  simulation  was  run,  with  results 
compared  to  observations  made  in  1999.  Inspection  of  simulated  and  observed 
contour  plots  revealed  that  simulated  TCE  concentrations  matched  well  to 
observed  data,  while  simulated  DCE  concentrations  were  underestimated  and 
simulated  VC  concentrations  overestimated.  The  value  of  E  for  TCE  (0.77) 
suggests  that  the  TCE  concentration  fit  is  good,  while  the  value  of  E  for  DCE 
(0.47)  indicated  that  fit  of  simulated  DCE  to  observed  values  was  relatively  poor 
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(Frankenberger  et  al.,  1999).  While  E  could  not  be  calculated  for  VC,  the  positive 
value  of  VC  bias  (B)  indicated  that  VC  was  being  over-estimated. 

While  chlorinated  ethene  advection  is  relatively  slow,  and  only  six  additional 
years  were  simulated,  we  found  that  BR3D  provided  good  predictions  of  TCE 
concentrations  and  fair  predictions  of  DCE  concentrations.  Due  to  a  lack  of 
quantifiable  data,  no  conclusions  were  drawn  regarding  VC  predictions.  From 
these  conclusions  regarding  model  performance,  we  made  the  inference  that 
TCE  and  DCE  fate  and  transport  could  be  simulated  by  an  advection/dispersion 
model  that  assumes  linear  equilibrium  sorption  and  first  order  decay. 

5.4  Areas  for  Further  Research 

1 .  In  Chapter  3,  sorption  was  assumed  to  be  a  linear  equilibrium  process.  To 
test  this  assumption,  the  Damkohler  number  was  calculated  for  the  site. 
The  Damkohler  number  I  (Dai),  defined  as  the  ratio  between  advection 
and  reaction  time  scales,  can  be  used  to  indicate  if  sorption  is  rate-limited 
or  in  equilibrium.  If  the  Dai  is  large,  then  advection  is  slow  in  comparison 
to  sorption,  suggesting  that  sorption  can  be  modeled  as  an  equilibrium 
process.  Conversely,  a  small  Dai  indicates  advection  is  fast  compared  to 
sorption,  and  that  sorption  would  be  best  modeled  as  a  rate-limited 
process.  The  value  of  Dai  depends  on  the  groundwater  velocity,  length 
scale  (typically,  diameter  of  the  porous  material),  and  the  sorption  rate 
constant.  Using  site  values  for  average  groundwater  velocity  and 
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literature  values  for  the  diameter  of  sand  (the  predominant  aquifer  material 
found  during  the  treatability  study)  and  a  sorption  rate  constant,  we 
calculated  Dai  =  1 .4x1  O'3.  Since  Da,  «1 ,  sorption  at  the  site  may  best  be 
modeled  as  a  rate-limited  process  instead  of  an  equilibrium  process. 

2.  While  it  is  hypothesized  that  DCE  and  VC  are  aerobically  oxidized 
downgradient  of  the  source,  the  lack  of  oxygen  data  for  the  site  precludes 
modeling  this  process.  Validation  of  a  site  with  more  comprehensive 
geochemical  data  might  instill  more  confidence  in  the  model’s  ability  to 
simulate  complex  biological  processes. 

3.  The  calibration  of  the  contaminant  transport  model  was  accomplished 
through  a  trial-and-error,  forward  process.  While  this  method  was 
adequate  for  the  purposes  of  this  research,  it  may  be  possible  to  obtain  a 
better  fit  to  the  observed  concentrations  using  a  nonlinear  simulation- 
regression  code,  such  as  UCODE  (Hill,  1998;  Gandhi  et  al.,  in  review),  to 
estimate  model  parameters. 

4.  This  research  modeled  the  natural  attenuation  of  chlorinated  ethenes  by 
reductive  dehalogenation.  As  additional  models  are  developed  to 
describe  the  processes  thought  to  be  taking  place  at  Site  LF-03,  the 
methodology  developed  in  this  thesis  can  be  used  to  measure  the 
predictive  ability  of  each  new  model.  The  resulting  best-fit  statistics  can 


74 


then  be  compared.  Hill  (1998)  provides  an  introduction  to  model 
comparison,  and  specifically  discusses  comparison  between  models  with 
different  numbers  of  input  parameters,  while  Legates  and  McCabe  (1999) 
discuss  the  necessary  methods  to  determine  the  statistical  significance  of 
the  coefficient  of  efficiency. 

5.  RMSE  and  E  are  overly  sensitive  to  extreme  values  due  to  the  squared 
differences  used  in  their  calculation.  Legates  and  McCabe  (1999)  discuss 
the  use  of  nonparametric  or  rank  correlation  methods,  as  well  as 
goodness-of-fit  statistics  that  rely  on  the  absolute  error  rather  than  the 
squared  error.  Such  methods  could  prove  useful  in  further  model 
validation  efforts. 
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